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Foreword 


A dissertation  prospectus  is  primarily  a research  proposal  and 
is  not  normally  published  or  released  outside  the  academic  institution 
involved.  This  prospectus  contains  an  unusual  amount  of  background 
material,  and  the  state  of  the  research  is  advanced  for  such  a docu- 
ment. This  prospectus  is  being  published  to  make  this  background 
material  and  the  analytical  methods  immediately  available. 

The  problem  of  gravity  modeling  for  terrestrial  inertial  naviga- 
tion is  discussed  at  length.  This  historical  and  technical  perspective 
provides  a survey  of  this  area  not  available  in  the  literature.  Since 
this  discussion  will  not  be  repeated  in  its  entirety  in  the  dissertation, 
publication  of  this  prospectus  will  provide  this  material  to  other 
interested  researchers. 

The  analytical  development  in  the  proposed  approach  (Part  III) 
includes  results  which  apply  to  more  general  model  evaluations  than 
the  subject  research.  These  methods  can  be  applied  to  current  studies 
of  geodetic  effects  on  navigation  accuracy.  Publishing  these  results 
will  make  them  immediately  available. 
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Abstract 

A historical  perspective  of  gravity  modeling  for  terrestrial  navi- 
gation is  presented.  The  traditional  ellipsoidal  model  is  explained, 
and  the  consequent  errors  are  discussed.  The  propagation  of  these 
errors  into  navigation  estimation  errors  is  presented.  A brief  survey 
of  advanced  modeling  methods  and  the  pertinent  theory  is  presented. 

The  system  design  problem  of  selecting  an  advanced  gravity 
model  is  presented  as  a scenario  to  motivate  the  proposed  research. 

To  address  this  problem,  a new  theoretical  analysis  .echnique  is 
developed.  This  technique  includes  the  effects  of  navigation  error 
propagation,  the  statistics  of  the  anomalous  field  (the  residual  after 
ellipsoidal  or  other  reference  field  modeling),  the  statistics  of  gravity 
survey  errors,  and  the  advanced  gravity  modeling.  These  effects  are 
combined  to  yield  a measure  of  system  performance  cost  as  reflected 
in  the  navigation  error  state  covariance  due  to  gravity  modeling  errors 
acting  alone. 

This  report  contains  82  references  to  the  open  literature  in  this 
subject  area. 


xiii 


I.  Purpose 


This  prospectus  defines  my  dissertation  research  topic  and 
provides  a basis  for  comments  or  approval.  To  this  end,  it: 

a.  Places  the  proposed  work  in  perspective; 

b.  Specifies  the  research  question  to  be  addressed; 

c.  Specifies  the  proposed  approach  and,  where  necessary, 
points  out  areas  which  will  not  be  addressed; 

d.  Specifies  the  method  for  confirming  the  final  result; 

e.  Outlines  the  Final  Report; 

f.  Describes  the  applicability  of  this  research  to  the  mission 
of  the  Air  Force  Avionics  Laboratory;  and 

g.  Discusses  the  aspects  of  this  research  perceived  to  be 
original. 
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II.  BackRTovind 


The  proposed  study  concerns  changes  to  the  current  methods  of 
accounting  for  gravitational  acceleration  in  inertial  navigation;  a brief 
discussion  of  the  state-of-the-art  is  in  order.  The  gravity  model*  is  a 
closed-form  mathematical  expression,  or  algorithm,  which  together 
with  a set  of  predetermined  parameters,  or  data,  define  a vector 
valued  function  which  approximates  the  Earth's  gravitational  accelera- 
tion throughout  some  three-dimensioned  region  of  operation.  Such  a 
model  is  an  integral  part  of  an  inertial  navigation  system  (INS).  Tech- 
nological advances  in  INS  instrument  designs  and  new  mission  require- 
ments call  for  a reappraisal  of  gravity  modeling  errors.  Improvements 
to  the  model  must  be  based  on  available  gravity  data;  on  our  knowledge 
of  gravitational  theory;  and  on  our  knowledge  of  the  time  spectral 
response  of  the  navigation  filter.  Of  course,  this  study  must  branch 
from  the  past  and  current  gravity  modeling  research.  We  begin  with 
an  assessment  of  where  gravity  modeling  evolution  has  led  to  date. 

The  term  "gravity  model"  can  be  m sleading.  The  term 
"gravitation"  will  always  refer  to  the  mass  attraction  (G)  effect  alone, 
whereas  "gravity"  will  usually  refer  to  the  combined  effect  of  mass 
attraction  and  Earth's  rotation  (£^). 


A.  The  Role  of  the  Gravity  Model  in  an  INS 

A gravity  model  is  a fundamental  component  of  any  precision 
INS.  We  shall  restrict  our  attention  to  "precision"  navigation.  An 
inertial  navigation  system  senses  inertial  translational  acceleration 
and  rotational  velocity,  tracks  attitude  by  integrating  the  rotational 
velocity,  and  calculates  translational  position  and  velocity  as  integrals 
of  acceleration  from  the  Newton's  Second  Law  equation:  Force  equals 
mass  times  acceleration.  The  mathematical  form  of  the  translational 
solution  algorithm  depends  on  both  the  coordinate  frame  in  which  the 
sensors  are  mounted,  instrument  or  platform  frame,  and  the  coordi- 
nate frame  in  which  the  computations  are  carried  out;  computational 
frame  (Ref  1).  Usually,  the  rotational  velocity  is  sensed  by  a set  of 
gyroscopes.  For  the  translational  acceleration,  sensors  called  accel- 
erometers are  used.  The  accelerometers  use  a test  mass  suspended, 
in  some  manner,  within  the  instrument  case  as  an  acceleration 
detector.  Thus,  the  accelerometer  only  senses  the  differential  accel- 
eration between  the  test  mass  and  the  instrument  case  which  is  mounted 
through  a suspension  system  to  the  body  of  the  vehicle.  This  differ- 
ential acceleration  is  called  specific  force,  f,  and  results  from  forces 
which  act  directly  on  the  vehicle  such  as  thrust  or  aerodynamic  lift. 
Gravitational  attraction  acts  simultaneously  on  the  teat  mass  and  the 

Examples  of  systems  not  included  are  those  which  do  not 
inertially  instrument  the  vertical  channel  (Reftl:109)  and  systems  with 
modeling  assumptions  which  would  mask  any  improvement  in  the 
gravity  model  (e.g.  Ref  2:133). 
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vehicle,  hence  it  cannot  be  sensed  by  the  accelerometers  (Ref  3s230- 
Z.'il,  and  Ref  1:2).  The  total  inertial  acceleration  is  the  sum  of  the 
sensed  specific  force  and  the  unsensed  gravitational  acceleration. 

Since  gravitational  attraction  is  not  sensed  directly,  an  INS  must  con- 
tain a gravity  model  to  complete  the  information  needed  for  the  trans- 
lational navigation  task. 

B.  Traditional  Modeling  Approaches 

The  nature  of  past  INS  gravity  models  for  airborne  use  (Ref  4: 
pp.  1-4,  6)  has  been  driven  by  two  predominant  factors.  First,  com- 
puter memory  space  and  computation  time  represent  a precious  system 
resource.  Other  mission  requirements  are  frequently  asynchronous  to, 
and  higher  priority  than,  the  navigation  function:  for  example  the 
attitude  control  function.  The  second  factor  is  the  existence  of  inertial 
instrument  errors.  The  sensed  specific  force  and  the  computed  gravity 
are,  symbolically  at  least,  summed  before  entering  the  first  level  of 
integration.  Thus,  there  is  little  incentive  for  creating  a gravity  model 
which  is  much  better  than  the  accelerometer  uncertainty  level.  These 
factors  have  led  to  the  use  of  simple  models  which  are  based  on  an 
intuitive,  but  effective,  idealization.  In  this  section  we  shall  discuss 
the  historical  basis  for  this  ideal  or  normal  gravity  model,  we  shall 
evaluate  the  effectiveness  of  this  ideal  field  as  an  approximation  to  the 
true  field,  and  we  shall  discuss  the  practical  approximations  to  this 
ideal  model  which  are  actually  used. 
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The  fundamental  shape  of  the  Earth  is  -well-described  by  an 
ellipsoid  (an  ellipse  of  revolution  about  its  minor  eixis  through  the 
poles).  This  notion  was  conceived  by  Newton  (Ref  8).  It  is  based  on 
the  concept  of  water,  the  seas,  not  flowing  on  an  equipotential  surface, 
mean  sea-level  for  example.  This  lack  of  flow  also  indicates  the  force, 
gravity  in  this  case,  is  perpendicular  to  the  surface  (hence  the  term 
"normal  field").  It  is  interesting  to  note  that  the  shape  of  a homo- 
geneous fluid  rotating  at  Earth  rate  would  be  an  ellipsoid  with  eccen- 
tricity very  near  that  of  the  reference  ellipsoids  used  in  geodetic 
surveys  and  systems  (Ref  9).  The  term  "ellipsoid  gravity  model" 
refers  to  a model  based  on  the  gravity  vector  being  perpendicular  to  the 
surface-approximating  ellipsoid. 

This  ellipsoidal  model  requires  only  three  parameters,  but  it  is 
impressive  in  terms  of  the  small  residual  error  between  model  gravity 

C 

and  true  gravity.  This  error  never  exceeds  400  mgal  (1  mgal  =•  lO'  m/ 
2 

sec  , approximately  1 Mg) --an  error  of  approximately  one  part  per 
2500.  Kayton  (Ref  5;Vol  II  p 289)  suggests  that  the  ellipsoidal  model  is 
adequate  for  systems  with  accelerometer  uncertainty  levels  greater 
than  20  mgal.  This  results  from  the  premise  that  we  view  acceler- 
ometer and  gravity  model  errors  in  a static  sense.  A more  precise 
measure  would  consider  the  time  spectral  properties  of  these  different 
system  error  sources.  We  must  understand  more  about  the  application 
of  the  ellipsoidal  model  before  this  point  can  be  pursued. 
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The  ellipsoidal  model  is  cumbersome  to  program  since  it  calls 
for  terms  not  necessary  for  any  other  navigation  function.  The  result 
is  approximations  of  the  approximation.  That  is,  computationally 
efficient  simplifications  are  used  in  place  of  the  complete  ellipsoidal 
formulae.  The  simplifications  are  developed  around  readily  available 
navigation  quantities,  and  this  set  varies  with  specific  navigation 
algorithm.  Two  examples  are  the  local-level  computational  frame  and 
the  inertial  computational  frame.  An  INS  which  computes  in  a local- 
level  (Ref  1:109)  coordinate  frame  tends  to  have  latitude  and  altitude 
available.  Hence,  approximations  are  made  based  on  these  inputs 
(Ref  4:1-11,  1-16,  1-23,  and  1-28).  For  Earth-centered  inertial*  (ait> 
known  as  geocentric  inertial,  ECI,  or  i-frame)  or  in  Earth-centered 
Earth-relative  (Geocentric  relative,  ECE,  or  e-frame),  the  "gravita- 
tional" field  resulting  from  an  ellipsoid  "gravity"  model  is  computed 
directly  from  the  Cartesian  components  of  the  position  vector.  For 
this  case,  the  formula  used  is  based  on  a truncated  series  expansion  of 
the  ellipsoidal  potential.  This  spherical  harmonic  model  requires  the 
definition  of  a number  of  coefficients  called  J -terms.  These  J-terms 
are  based  on  satellite  and  ground  survey  gravity  data.  This  approach 
is  well-documented  (see  Ref  1:49-53;  Ref  5:309-313;  Ref  6:195-225; 

Ref  9:Appendix;  Ref  10;App.B;  Ref  IhApp.E;  Ref  12:35-43;  Ref  13:36- 
42;  Ref  14:20-35;  Ref  15:139-142;  and  Ref  16:419-423).  We  shall  return 

Typically  these  frames  are  centered  at  the  Earth's  center  of 
mass  (Ref  5:Vol  I,  pp  6-9  and  Vol  D.,  pp  288-289);  the  i-frame  is  non- 
rotating with  respect  to  the  distant  stars. 
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to  this  particular  model  later,  but  for  now,  consider  the  results  of 
these  approximations  to  the  ellipsoid  gravity  model. 

TLo  standard  for  comparison  of  modeling  simplifications  is 
the  WGS  72  ellipsoid  (Ref  17).  The  local-level  models  vary  as  much 
as  20  mgal  from  this  reference.  This  result  is  due  to  the  altitude 
correction  simplifications.  The  spherical  harmonic  model  can  be  more 
precise;  it  reproduces  the  ellipsoidal  model  to  within  0.  1 mgal  in  one 
example  (Ref  4:1-5). 

Direct  comparison  of  these  simplified  ellipsoid  models  with 
true  gravity  have  not  been  made.  Geodetic  gravity  surveys  measure 
gravity  on  the  Earth's  surface  and  extrapolate  to  a reference  ellipsoid. 
This  measure  of  gravity  at  the  ellipsoid  surface  can  be  compared 
directly  to  the  ellipsoid  gravity  model  value  (Refs  9,  14,  and  17).  A 
good  ellipsoid  can  result  in  global  maximum  error  of  approximately  400 
mgal.  The  ellipsoid  can  also  be  varied  to  produce  a better  local  fit  in 
some  region  of  interest.  The  nature  of  these  errors  in  modeling  is 
important  and  merits  further  discussion. 

The  Nature  of  Modeling  Errors 

The  Earth's  gravity  field  varies  with  position  and  with  time. 

The  time  variations  can  be  subdivided  into  categories: 

a.  The  Earth's  rotation  with  respect  to  inertial  space, 

b.  Celestial  bodies  (primarily  the  Sun  and  Moon),  and 

c.  Earth's  mass  redistributions  (primarily  tides). 
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The  effects  due  to  b.  and  c.  are  on  the  order  of  tenths  of  a mllli-gal 
(Ref  18:76-77).  The  effects  due  to  a.  are  on  the  order  of  hundreds  (not 
hundredths)  of  mllli-gals.  Although  modeling  the  effects  due  to  b.  and 
c.  would  be  relatively  simple,  this  effort  would  be  pointless  without  a 
major  resolution  of  the  a.  effects.  So,  we  shall  restrict  our  attention 
to  the  Earth  rotation  effect  which  is  simpler  from  an  Earth-relative 
observer's  point  of  view.  Since  we  are  ignoring  the  tidal  and  celestial 
body  effects,  an  Earth-relative  observer  would  see  a static  gravity 
field  at  each  point.  The  true  (time-average.  Earth-relative)  gravity 
field  is  then  a vector  function  of  an  Earth- relative  position  vector 
(e-frame  for  convenience).  Symbolically,  we  are  stating  that 

G®(r,  t) » G®(r®)  ( 1) 

The  model  we  seek  is  a closed-form  mathematical  relationship  either 
Gp^(r^)  to  approximate  gravitation  ^(£®)  or  to  approximate 

gravity  £(r®). 

Going  back  to  our  inertial  observer,  we  may  mathematically 
express  what  he  will  see  in  terms  of  the  above  function  definitions. 

G'(r®)  = G®(r®)  = c'gG®(C?r')  (2) 

where  is  the  e-frame  to  i-frame  coordinate  transformation  matrix; 
Ci  is  its  inverse.  Elements  of  are  trignometric  functions  of  the 
product  of  time  and  Earth  rotation  rate,  w^^t  (Ref  1:36  for  example). 

So,  for  a fixed  r^,  the  gravity  field  is  periodic  at  Earth  rate.  Then, 
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1 0 
C^{’)  is  an  explicit  function  of  time  through  and  an  implicit  function 

of  time  through  £®(t).  With  our  assumptions,  we  consider  G®(*  ) to  be 
only  an  implicit  function  of  time  via  r_®(t). 

We  may  now  formally  define  the  gravitation,  or  gravity,  model- 
ing errors  by 


tv 

(3) 

Since, 

& = G X (wieX  r). 

(4) 

we  also  have 

-iiieXiiLie  X£  - [G^(r®)  - Wie  XoLig  X£ ]. 

Hence,  5£(£®)  * - G(r®).  (5) 

The  gravitational  and  the  gravity  modeling  errors  are  the  same.  This 
process  is  depicted  in  flow-chart  form  in  Figure  1. 

D.  Propagation  of  Errors  in  an  INS 

The  resulting  error  vector,  5^,  is  an  explicit  time  function  in 
the  i-frame  and  an  implicit  function  of  time  in  the  e-frame.  Since  • ) 
models  the  Earth- relative  time-average  gravitational  field,  ^®(*),we 
must  view  5£.®(*)  as  the  time  average  gravitational  error  for  each 
argument  Having  defined  this  error,  we  now  turn  our  attention  to 

the  effect  it  has  on  navigation  performance.  We  need  .o  understand  how 

* 

these  errors  enter  the  navigation  algorithm,  what  in  mathematical 


terms  the  algorithm  does  to  these  errors,  and  the  nature  and  extent  of 
corruption  of  navigation  estimates. 

The  INS  solves  for  position  and  velocity  by  a double  integration 
of  an  estimate  of  acceleration.  This  estimate  is  formed  by  summing 
the  specific  force  sensed  by  accelerometers  with  the  model-computed 
gravitational  acceleration.  We  must  clearly  distinguish  between  the 
modeling  error,  3£,  and  the  INS  propagation  error  which  results  from 
an  error  in  the  INS  position  estimate.  This  propagation  error  can  be 
formally  defined  as  the  gravitation  at  the  estimated  position,  ^(r),  less 
the  gravitation  at  the  true  position,  G(r).  With  this  distinction  between 
the  different  types  of  gravitational  modeling  errors,  we  may  proceed 
with  an  INS  error  analysis. 

We  are  interested  in  INS  estimates,  so  we  need  to  define  posi- 
tion and  velocity  error  quantities.  Let 


« r(t)  = r(t)  - r,t) 

(6) 

= £{t)  - r(t)  . 

(7) 

The  exact  form  of  the  relationships  between  ir_,  and  6jr  depends  on 
the  specific  navigation  mechanization  and  computational  approach. 

Such  matters  as  independently-sensed  velocity  feedback  and  vertical 
channel  altimeter  aiding  introduce  analysis  complexities  that  would 
obscure  the  fundamental  point  we  are  pursuing  at  this  time.  These 
issues  must  be  dealt  with,  but  we  shall  table  them  for  now.  We  can  see 
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the  basic  interrelationship  by  considering  a first-order  perturbation  of 
the  navigation  differential  equation 

i = l + G(r®),  (8) 

whereof  is  the  true  specific  force. 

The  specific  form  of  (8)  depends  on  the  computational  coordinate  frame. 
Reference  1 presents  a full  description  of  this  topic  and  introduces  the 
notation  used  herein.  References  2 and  19  are  extensions,  corrections 
and  clarifications  of  this  work.  Many  error  sources  must  be  included 
to  perform  a complete  navigation  error  analysis.  We  are  restricting 
our  scope  to  gravity  modeling  errors  and  their  propagation  (e.  g. 
if_  =^).  The  first-order  differential  equation  from  (8)  becomes 

if  = «^(r®)  + r (r®)  « r (9) 

where  r(£®)  is  the  gravitation  tenr.>r  (see  Appendix  A). 

This  equation  describes  the  propagation  of  both  types  of  gravity 
errors:  gravity  modeling  errors  given  by  j£{_r®)  and  the  error  in 
calculating  gravity  with  an  error  in  the  position  estimate.  Since  we 
assume  no  velocity  or  altitude  aiding,  we  may  obtain  by  a double 
integration  of  (9).  Appendix  B shows  the  inherent  instability  of  the 
"vertical  channel"  of  (9).  In  practice  altimeter  aiding  is  used  to 
stabilize  this  channel  and  this  will  complicate  the  analysis.  It  will  suit 
our  purposes  to  consider  only  the  horizontal  channels  for  a simple 
example  of  error  propagation.  Again,  Appendix  B develops  the 
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homogeneous  portion  of  (9)  for  one  horizontal  coordinate.  Adding  the 
modeling  error  back  in  yields 

..  2 

6x  = 5g„  - (u  ix  (10) 

s 

where  oig  is  the  Schuler  rate  given  by  yJT/^ • 

This  undamped,  second-order  error  differential  equation  has 
well-known  characteristics.  Sinusoidal  steady-state  analysis  shows  a 
high  sensitivity  to  time -frequencies  near  this  Schuler  rate.  The  steady- 
state  analysis  fails  for  inputs  of  at  the  Schuler  rate  since  the  par- 
ticular solution  to  (10)  is  unbounded.  This  case  can  be  anal^ised  in  the 
time  domain  and  can  give  some  insight  into  the  problemb  that  gravity 
modeling  errors  might  cause.  Reference  20  provides  data  on  the  vari- 
ation of  gravity,  or  gravitation,  from  a reference  ellipsoidal  field 
along  a transcontinental  path  across  the  35°  latitude  line  of  the  United 
States.  The  horizontal  component  is  recorded  as  the  angle  true 
gravity  is  deflected  from  the  ellipsoidal  gravity  in  East-West,  prime, 
and  North-South,  meridian,  directions.  The  root- mean  -square  value 
for  each  direction  is  near  five  arcseconds  (Ref  21),  which  translates 
to  approximately  25  mgal.  Suppose  the  driving  term  in  (10)  is 
sinusoidal  in  its  argument  for  horizontal  motions.  Let  d be  the  spatial 
period  so 

This  example  is  contrived  to  show  the  worst  accuracy  degrada- 
tion, however  it  provides  a simple  closed-form  demonstration  of  some 
basic  concepts. 


(11) 


Now  suppose  the  vehicle  travels  at  constant  velocity  such  that 
Vjj.  = ojigd/2Jr.  Express  x as;  x = v^^t  + 0 = (o)gd/27r)t.  Then 
5gx  = 25  sin  (o)gt)  mgal. 

Identify  ix  as  and  assume  zero  initial  conditions  on  (10). 

We  can  solve  the  particular  problem  now  by  LaPlace  transforms  and 
find  ^^5t(t)  = (25  mgal)  t sin  (Wgt).  (13) 

Note  that  the  amplitude  of  the  sinusoid  grows  linearly  with 
time--an  unbounded  response.  Take  for  example  the  air  launch  of  an 
ICBM  as  a way  of  gauging  our  concern  with  (13).  We  know  (Ref  16:305) 
that  at  ICBM  ranges  a one  foot-per-second  velocity  error  can  cause  a 
one  nautical  mile  miss.  At  this  rate,  we  would  only  have  two  minutes 
of  cruise  time  before  (13)  would  indicate  error  divergence  with  the 
potential  for  a 0.  1 n.  mi.  miss  from  gravitational  modeling  effects  alone. 
The  point  of  this  much-contrived  example  is  that  there  are  conceivable 
applications  where  present  gravity  models  are  inadequate.  Also,  it 
should  be  clear  that  the  extent  to  which  our  estimates  are  corrupted  by 
this  gravity  noise  depends  on  the  spatial  frequency  of  the  gravity  model- 
ing errors,  on  the  time  frequency  response  of  the  INS,  and  on  the 
mission  velocity  which  translates  the  spatial  gravity  modeling  error 
function  into  the  time  domain. 

E.  Magnitude  of  Resulting  INS  Errors 

With  this  background,  let  us  consider  more  complete  studies  of 
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gravity  modeling  effects  on  INS  estimates.  The  open  literature  con- 
tains two  distinct  types  of  studies  on  navigation  accuracy  with  gravity 
modeling  errors:  statistical  and  deterministic.  The  statistical  studies 
(Refs  21,  22  and  23)  are  covariance  analyses  or  Monte  Carlo  analyses 
based  on  a statistical  model  of  the  gravity  modeling  disturbance.  The 
deterministic  studies  (Refs  24  zind  25)  are  simulation  analyses  which 
include  detailed  local  gravity  models  for  the  region  of  interest.  The 
statistical  studies  provide  parametrically  the  expected  accuracy  degrada 
tion,  whereas  the  deterministic  studies  give  specific  case  data  which 
include  error  time  histories. 

1,  Statistical  Analyses.  Historically  the  statistical  studies 
appear  first.  The  gravity  modeling  error  is  itself  statistically  modeled. 
The  disturbance  vector,  6^^,  is  not  necessarily  modeled  directly;  in 
Reference  21,  the  deflection  angles,  | and  t),  are  modeled  as  the  outputs 
of  spatial  linear  systems  driven  by  white  gaussian  noise  (Ref  26), 
Statistical  model  details  vary,  and  a great  deal  of  research  has  been 
devoted  to  defining  a statistical  model  which  is  consistent  with  both 
observed  gravity  data  and  with  gravitational  field  theory.  The  statisti- 
cal model  issues  are  discussed  in  Appendix  C;  our  concern  here  is  the 
effect  of  gravity  modeling  errors  on  navigation  accuracy. 

The  Levine  and  Gelb  paper- (Ref  21)  is  the  classic  in  the  statisti- 
cal studies  category.  Their  approach  is  a steady-state-covariance 
analysis  (Ref  27)  which  requires  the  total  problem  to  be  cast  in  the 
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•form  of  a linear,  time-invariant  stochastic  differential  equation.  The 
statistical  model  for  the  gravity  error  is  based  on  an  exponential  auto- 
correlation function,  for  example 

(X)  = e-  (U) 

where  | is  the  meridian  vertical  deflection, 

X is  the  shift  distance, 
is  the  variance  of  %,  eind 

d|  is  the  correlation  distance. 

2 

With  (14),  (T^  and  dg  completely  specify  the  statistical  model.  Gravity 
disturbance  components  are  assumed  uncorrelated,  so  | and  q (the 
prime  deflection)  are  uncorrelated.  This  assumption  statistically 
uncouples  the  horizontal  channels.  The  INS  model,  also,  dynamically 
uncouples  these  channels.  Rate  feedback  from  a non-inertial  sensor  is 
assumed  and  the  error  propagation  is  governed  by  (10)  with  the  addition 
of  the  rate  feedback  term: 

6x  = -K  6x  - Wg  6x  - 6gjj  ’ (15) 

where  K is  the  feedback  gain  for  the  non-inertial  velocity  damping. 
Recall  that  the  disturbance  term  is  a function  of  Earth-relative  position. 
To  implement  (15)  we  must  convert  the  disturbance  term  to  the  time 
domain.  This  is  accomplished  through  two  steps:  model  (14)  as  a 
first-order,  linear,  time-invariant  stochastic  differential  equation  in 


i 
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the  argument  x and  use  velocity  to  convert  x,  hence  (14),  to  the  time 
domain.  First, 


i-lw+q,  (X) 

dx  5 


(16) 


where  q>(x)  is  a zero-mean,  white,  gaussian  noise  process  with  auto- 
^ 20-^ 

correlation  function  of  L S(x).  5(x)  is  the  Dirac  delta  function  which 

satisfies  I _ Q for  all  X 0, 


f 6(x) 

-€ 


dx  = 1 for  all  € > 0 . 


(17) 


The  entire  analysis  is  performed  in  a quasi-static  sense  with  a constant 
velocity  assumed  then  varied  parametrically.  Then, 


and 


dx  = Vjj  dt. 


(18) 


6[x(t)]  = 6 (vx  t)  = ~L  6 (t) 


(19) 


So, 


(20) 


where  q^(t)  is  zero-mean  white  gaussian  noise  with  autocorrelation  of 


2 a! 


(?) 


6{t),  Because  deflection  angles  are  small  we  assume 


66x  = U- 


(21) 


With  (21),  equations  (15)  and  (20)  form  the  required  stochastic  time 
differential  relationship  for  the  steady-state  covariance  analysis.  In 
retrospect,  it  is  interesting  to  compare  the  transition  from  spatial,  (11), 
to  time,  (12),  domain  in  the  simple  closed-form  example  above  to  the 
transition  in  this  stochastic  analysis  from  (16)  to  (20). 

With  this  system  of  differential  equations  driven  by  white 
gaussian  noise,  the  steady  state  covariance  of  navigation  position  and 
velocity  errors  can  be  found  using  LaPlace  transform  techniques  on  the 
system  covariance  matrix  equation.  The  structure  of  the  Levine  and 
Gelb  analysis  is  easier  to  summarize  than  the  results.  Velocity,  corre- 
lation disteince,  and  INS  damping,  K,  are  varied  parametrically.  The 
independent  variations  of  v^  and  d^  are  tinnecessary  since  in  the  final 
results  they  always  appear  in  the  ratio 

Pg  =Vx/d^  (22) 

This  result  is  important  in  its  own  right;  this  implies  that  the  same 
rms  navigation  error  results  from  an  increase  in  velocity  as  in  a 
decrease  in  correlation  disteince.  In  the  limit,  increased  velocity  or 
decreased  correlation  distance  cause  the  gravity  disturbance  to  approach 
white  noise  (which  is  precisely  what  zero  correlation  distance  would 
mean  in  the  spatial  domain)  with  the  exponential  function  approaching 
the  delta  function  and  correlation  time  (1/P^)  approaching  zero.  In  the 
limit  in  this  direction,  both  position  and  velocity  rms  errors  approach 
zero. 
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For  the  other  extreme,  consider  either  velocity  decreasing  or 
correlation  distance  increasing  (|3^  approaching  0).  Then,  the  disturb- 
ance looks  more  like  a constant.  Since  we  have  a non -inertial,  and  for 
now  presumably  perfect,  velocity  sensor,  the  steady-state  velocity 
error  goes  to  zero.  The  position  error  approaches  the  constant  value 
that  reflects  the  position  offset,  incurred  while  nulling  the  velocity, 
required  to  exactly  cancel  the  gravity  disturbance.  From  (15),  this 
offset  can  be  seen  as  the  gravity  disturbance  divided  by  the  square  of 
the  Schuler  rate. 

While  the  position  rms  error  is  mzucimum  for  this  near-constant 
disturbance,  the  velocity  rms  error  behaves  quite  differently.  As  dis- 
cussed, the  velocity  rms  error  approaches  zero  for  the  cases  of  |3^ 
approaching  either  zero  or  infinity.  This  result  is  predicated  on  o-^ 
staying  finite  and  some  further  analysis  is  really  required  if  we  want  to 
let  8gx(t)  approach  white  noise.  The  velocity  rms  error  per  unit  of 
standard  deviation  has  a maximum  when  equals  the  Schuler  rate  ojg. 
The  8gjj.  spectral  characteristics  are  best  aligned  with  the  navigation 
filter  under  these  velocity/correlation  distance  conditions.  Hence,  in 
studying  gravity  error  propagation,  we  need  to  pay  attention  to  this  type 
alignment  since  the  error  sensitivity  is  highest  there.  In  particular,  if 
we  decompose  the  gravity  disturbance  into  a Fourier  series  from  the 
spatial  domain,  we  get  a sum  of  sinusoidal  terms  similar  to  (11).  The 
velocity  interpolation  can  then  be  used,  as  in  (12),  to  convert  to  the 
time  domain.  Then,  we  can  identify  those  frequencies  which  will  cause 


19 


the  greatest  navigation  rms  errors. 

Other  observations  can  be  made  on  the  various  interplays 
between  the  statistical  model  and  the  INS  dynamics;  however  these 
results  are  somewhat  intangible.  What  we  need  is  a perspective  on 
how  realistic  gravity  disturbances  compare  with  other  navigation  error 
sources.  The  traditional  way  of  presenting  these  data  is  an  error 
budget  which^for  a specified  set  of  hardware  and  a specified  mission^ 
allocates  the  expected  error  into  compartments  labeled  by  the  various 
recognized  error  contributors.  This  more  satisfying  approach  was 
taken  by  Nash,  D'Appolito,  and  Roy  in  Reference  23. 

For  a given  polar  flight  profile,  a gravity-error-alone  analysis 
shows,  over  the  first  phase  of  the  mission,  a 0.03  n.  m,  position  error 
per  hour  of  flight  per  arcsecond  rms  deflection  of  the  vertical.  Next, 
an  integrated  system  based  on  a 0.  5 n.  m.  /hr.  INS  hardware.  Doppler 
radar  velocity  aiding  and  Loran  position  aiding  were  included  in  the 
integrated  system  analysis.  A severe  gravity  disturbance  model  was 
used  with  a 20  arcsecond  rms  deflection  and  20  n,  m.  correlation 
distance.  In  the  final  analysis,  the  errors  due  to  vertical  deflections 
were  substantial.  Approximately  seventy-five  percent  (75%)  of  the 
radial  velocity  error,  thirty  percent  (30%)  of  the  radial  position  error, 
and  fifteen  percent  of  the  heading  error  were  p*,i.ributed  to  the  vertical 
deflections.  Even  if  the  input  deflections  were  reduced  to  a more 
globally  representative  rms  level,  say  seven  aircseconds,  the  resulting 
navigation  errors  would  constitute  a significant  and  irreducible  error 
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term  with  the  traditional  gravity  model.  More  significantly,  the  verti- 
cal deflection  induced  errors  would  be  on  the  same  level  as  the  gyro  and 
accelerometer  induced  errors.  A significant  improvement  in  the 
inertial  sensors,  say  to  a 0. 1 n.  m.  /hr  system,  would  not  correspond- 
ingly improve  integrated  system  accuracy  since  the  gravity  modeling 
error  would  remain. 

2.  Deterministic  Analyses.  From  both  parametric  and  specific 
mission  studies,  the  statistical  approaches  have  shown  the  nature  and 
extent  of  navigation  errors  induced  by  gravity  modeling  errors.  The 
deterministic  approach  allows  a means  of  comparing  these  statistical 
results  with  a gravity  field  which  is  very  nearly  a duplication  of  the 
actual  field  in  some  area  of  the  world.  Also,  by  forming  a truth  model 
to  simulate  the  actual  field  we  can  simulate  a navigation  mission  and 
produce  error  time  histories.  These  data  complement  the  average,  or 
expected,  data  from  statistical  analyses. 

Chatfield,  Bennett,  and  Chen  presented  such  an  analysis  in 
Reference  24.  For  a specific  flight  path  in  the  western  United  States, 
they  constructed  an  extensive  point- mass  anamolous  gravity  model  from 
available  measured  gravity  data.  We  shall  discuss  this  modeling  tech- 
nique later;  for  now,  it  is  sufficient  to  note  that  this  point-mass  model 
together  with  the  reference  ellipsoid  model  form  a much  closer  approx- 
imation of  the  true  field  in  a deep  volume  of  space  surrounding  the 
flight  path.  The  comparison  of  model  gravity  with  the  surface  gravity 
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data  shows  residuals  within  a mgal  directly  over  the  grid  midpoints. 

With  this  near-perfect  model,  simulated  aircraft  test  runs  were 
made  assuming  both  aided  and  unaided  INS.  The  point  mass  truth  model 
was  used  to  generate  true  position  and  velocity.  The  simulated  IN»S 
used  only  the  ellipsoidal  model.  Differences  between  the  two  gravity 
calculations  were  as  high  as  50  mgals  during  the  simulated  missions. 
The  differences  in  position  and  velocity  between  INS  estimates  and  the 
truth  model  output  represent  the  INS  estimation  errors.  These  errors 
are  indicative  of  what  one  could  expect  in  an  actual  flight. 

The  unaided  INS  errors  reached  peaks  of  near  2000  ft  in  posi- 
tion, 2 ft/sec  in  velocity,  and  15  arcseconds  in  heading.  The  aided  INS 
case  assumed  position  checkpoints  periodically  spaced  along  the  flight 
path  and  assumed  doppler  radar  velocity  measure.  Since  only  gravity 
anomaly  errors  were  simulated,  these  aids  were  noise-free  and  sub- 
stantially improved  performance.  The  position  errors  were  suppressed 
below  400  ft,  while  velocity  error  stayed  under  1 ft/sec.  The  heading 
angle,  while  smaller  than  in  the  unaided  case  most  of  the  time,  did 
reach  a peak  value  near  20  arcseconds.  A cursory  comparison  of  the 
error  time  histories  shows  for  the  unaided  case  that  the  statistical 
methods  of  Levine  and  Gelb  (Ref  21)  would  have  accurately  character- 
ized the  rms  errors.  The  time  histories  of  errors,  of  course,  would 
not  be  available  from  a purely  statistical  approach. 

Reference  25  contains  results  from  a continuation  of  the  Refer- 
ence 24  studies.  These  results  are  interesting  because  they  point  out 
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a problem  associated  with  gravity  modeling  beyond  the  INS  estimate 
errors.  In  a flight  test  of  an  INS,  one  goal  is  to  assess  INS  accuracy 
and  to  also  determine  the  accuracy  of  individual  components  and  sub- 
systems. The  usual  approach  is  to  structure  the  analysis  along  filter- 
ing theory  lines.  Component  error  sources  are  modeled  and  integrated 
with  the  INS  error  dynamics  model  to  form  a system  model.  Check- 
point and  other  tracking  data  give  a measure  of  the  INS  error  at  various 
points  during  the  mission.  Post  flight  data  processing  typically  employs 
regression  analyses  to  find  the  component  error  model  parameters 
which  best  fit  the  observed  data  in  the  sense  that  the  error  residuals  are 
minimized  in  some  way.  These  error  parameters  can  only  be  identified 
by  the  spectral  properties  of  their  resulting  INS  errors.  When  gravity 
model  errors  are  not  included  in  the  above  analysis,  the  gravity  model- 
ing induced  errors  spill  over  spectrally  into  this  component  parameter 
identification  process  thus  corrupting  the  resuPs.  This  problem  has 
been  ignored  heretofore  since  INS  component  errors  were  large  in  com- 
parison. This  study  (Ref  25)  was  conducted  for  Holloman  AFB  INS  test- 
ing and  is  evidence  that  the  day  has  arrived  when  we  can  no  longer 
ignore  this  factor.  Plans  for  testing  state-of-the-art  inertial  systems 
are  beginning  to  include  the  requirement  for  detailed  local  gravity  data 
for  use  in  the  regression  analyses. 

F,  Impetus  for  Improvements 

We  have  seen  the  nature  of  gravity  modeling  errors  and  the 
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magnitude  of  INS  errors  they  cause.  It  is  easy  to  understand  why  a 
detailed  new  model  might  be  used  in  a flight  test  environment  where  data 
purity  is  important.  One  could  question  the  need  for  model  refinements 
since  the  errors  induced  are  not  unacceptable  for  most  navigation  appli- 
cations. The  impetus  comes  from  the  potential  military  application- - 
where  the  impetus  for  inertial  navigation  originated.  In  the  delivery  of 
weapons,  any  error  diminishes  weapon  system's  effectiveness,  so  these 
gravity  induced  errors  cannot  be  ignored.  The  self-contained  nature  of 
inertial  navigation  virtually  assures  a continued  dependence  even  with 
advanced  radiometric  navigation  systems  available.  The  evolution  of 
strategic  force  concepts  provides  the  prime  motivation  for  the  increased 
emphasis  on  refining  the  traditional  gravity  modeling  techniques. 

The  United  States'  lead  in  quality  of  ICBMs  is  expected  to 
decline  (Ref  26)  in  the  next  decade.  One  proposed  response  to  a growing 
counterforce  threat  is  to  mobilize  the  currently  silo-based  ICBMs  (Refs 
27  and  28).  Leaving  the  static  silo  environment  for  a dynamic,  air 
launch,  or  intermittently  dynamic,  revetment  launch,  environment  will 
certainly  decrease  the  expected  accuracy  of  this  force  (Ref  29).  One 
cause  for  this  loss,  is  the  fact  that  silo-based  missiles  are  targeted 
with  the  aid  of  extensive  point-mass  models  for  the  launch  region  (Ref 
30).  While  this  loss  could  eventually  be  recovered  by  modeling  all  the 
possible  mobile  launch  regions,  we  still  have  to  deal  with  the  problem 
of  initial  conditions  for  the  missile  navigation  algorithm. 


24 


The  silo-based  missile  INS  can  be  initialized  with  the  known 


position  and  with  an  (expected)  zero  Earth-relative  velocity  prior  to 
launch  activities.  In  any  dynamic  launch,  these  data  must  be  provided 
by  navigation.  Our  earlier  example  of  the  nature  of  INS  error  propaga- 
tion indicated  that  in  that  contrived  case  only  minutes  would  elapse 
before  gravity  modeling  errors  would  blunt  the  weapon  system's  effec- 
tiveness. The  deterministic  study  (Ref  24)  of  the  last  section  indicates 
that  for  that  case,  even  with  many  external  aids,  the  velocity  errors 
can  grow  to  the  one-foot-per-second  level.  At  ICBM  ranges  this  error 
has  the  potential  to  cause  a 1 n mi  miss  (Ref  16:305).  From  Reference 
26,  this  type  of  accuracy  virtually  eliminates  such  a weapon  system 
from  use  against  targets  hardened  to  nuclear  attack. 

The  ICBM  concern  is  not  the  only  one.  Proposed  air-launched 
cruise  missiles  (Ref  31)  will  encounter  similar  problems.  The  potential 
degradation  may  be  greater  due  to  the  longer  navigating  flight  time  of 
the  missile.  Also,  one  must  consider  the  counterforce  mission  these 
cruise  missiles  are  postulated  to  perform. 

From  these  uniquely  military  applications,  we  can  see  a need 
to  refine  gravity  modeling.  The  nature  of  the  refinements  must  be 
predicated  on  the  real-time  and  on-board  need  for  the  data.  The 
dynamic  nature  of  future  missions  is  likely  to  preclude  the  extensive 
ground-based  pre-mission  targeting  which  compensates  for  these 
gravity  effects  today.  So,  whatever  new  methods  evolve,  they  must  be 
computationally  efficient;  yielding  the  greatest  accuracy  improvement 
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for  the  investment  of  on-board  computer  storage  space  and  execution 
time. 


G.  Improvement  Approaches 

There  are  many  current  and  past  activities  that  either  directly 
or  indirectly  suggest  approaches  to  this  problem.  The  "software" 
approaches  fall  into  several  categories.  The  statistical  approach 
simply  employs  filter  theory  to  fcrm  an  estimate  of  the  gravity  distur- 
bance. The  point- mass  model  is  an  example  of  finite -element 
approaches  which,  in  essence,  build  mass  distribution  perturbat.on 
models.  Another  software  approach  simply  employs  fundamental 
potential  integral  relationships.  These  integrals  are  approximated 
using  observed  gravity  data.  Then  transform  techniques  can  be 
employed  to  decrease  the  computational  burden.  Aside  from  these 
approaches,  new  inertial  instruments  are  being  developed  for  gravity 
mapping  (Ref  32).  As  McKinley  has  noted  (Ref  29),  it  is  ironical  that 
the  instrument  (gradiometer)  developed  for  the  mapping  application  can 
simultaneously  ma’^e  the  mapping  unnecessary.  We  need  to  delve  into 
all  these  areas  before  addressing  a plan  of  attack  on  the  model  refine- 
ment problem. 

1.  Statistical 

Statistical  methods  can  be  employed  pre-mission  (a  priori) 
or  during  the  mission  (real-time).  Pre-mission  data  preparation  might 
include  forming  a mission-dependent  ellipsoid  model  to  lower  the  rms 
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anomalous  gravity  for  that  specific  mission  trajectory  or  area. 

Another  a priori  approach  is  to  define  the  navigation  mission  in  great 
detail,  simulate  the  trajectory  with  an  accurate  gravity  model,  and  off- 
set the  initial  conditions  of  the  navigation  filter  to  compensate  for  the 
expected  errors.  Such  a priori  methods  are  open-loop  in  that  the  actual 
trajectory  may  vary  considerably  from  the  planned  one  due  to  environ- 
mental effects  or  other  modeling  errors.  Changing  the  INS  algorithm 
to  reject  the  short  spatial  wavelength  gravity  noise  is  unacceptable 
because  we  would  also  reject  any  measured  acceleration  which  fell  into 
the  same  time  spectrum.  The  real-time  role  for  statistics  has  little 
more  potential. 

In  real-time,  we  can  use  statistical  methods  to  predict  the 
near-term  influence  of  anomalous  gravity  on  navigation  estimates.  To 
implement  such  a scheme,  we  need  external  navigation  fixes  which  allow 
us  to  observe  the  navigation  errors  periodically.  With  anomalous  grav- 
ity modeled  as  a Gauss-Markov  process  (see  Appendix  C)  in  our  Kalman 
filter,  the  estimated  states  associated  with  this  model  will  allow  us  to 
statistically  predict  our  near-term  anomalous  gravity  effects.  This 
approach  has  been  considered  for  INS  flight  tests  (Ref  25)  as  a means 
of  removing  the  effects  of  gravity  noise  from  INS  component  perform- 
ance estimates. 

These  statistical  methods  have  limited  scope  because  they  do 
not  take  advantage  of  all  available  gravity  data  except  in  an  average 
sense  from  the  empirical  autocorrelation  function  (see  Appendix  C). 
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The  following  methods,  in  one  form  or  another,  use  the  gravity  data. 


2*  Finite  Element  Methods 

The  term  "finite  element"  to  some  connotes  a recursive 
approximation  based  on  a tri-diagonal  Grammian  matrix.  We  use  the 
term  to  denote  all  gravitation  or  geopotential  modeling  which  are  based 
either  on  a finite  partition  of  the  geoid  surface  for  integration  approxi- 
mation, on  a finite  number  of  mass  distribution  elements,  or  on  a finite 
set  of  local  approximating  (interpolating)  f\inctions.  Of  all  the  methods 
we  shall  discuss,  only  the  point-mass  method  has  ever  been  used  for 
INS  aiding--and  this  aiding  was  pre-mission^not  in-flight  (Ref  30),  So, 
these  methods  are  discussed  because  they  have  the  capability  to  form 
improved  INS  gravity  models.  All  of  these  models  have  been  proposed 
to  support  the  basic  task  of  Physical  Geodesy:  defining  the  geoid  (Refs 
8 and  34). 

Since  we  already  have  Legendre  polynomial  spherical  harmonic 
functions  as  a spanning  set  of  fvinctions,  one  might  ask  why  we  do  not 
fill  out  the  coefficients  to  the  degree  and  order  necessary.  The  trun- 
cated spherical  harmonic  series  forms  a finite  element  model  by  the 
above  definition,  so  we  shall  treat  it  as  such  and  discuss  its  merits  as 
an  improvement  candidate.  Since  this  model  is  global  in  scope,  it 
requires  global  data  for  complete  coefficient  identification.  The  Earth 
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Local  models  can,  however,  be  created  by  using  restricted- 


area  data. 


is  large  and  we  must  have  informatifT  in  our  model  of  relatively  short 
spatial  wavelength.  Shannon’s  sampling  theorem  tells  us  to  collect  data 
on  a grid  finer  than  the  shortest  wavelength  that  we  wish  to  represent-- 
such  a survey  is  not  economically  or  politically  feasible.  With  extensive 
new  data  from  GRAVSAT /GEOPAUSE,  Koch  (Ref  35)  points  out  the 
potential  for  rather  severe  aliasing  at  the  order-and-degree  12  trvinca- 
tion  level.  The  message  is  clear,  if  you  want  to  significantly  improve 
the  gravity  model,  you  should  concentrate  your  model  and  your  survey 
in  the  local  area  of  operation. 

As  mentioned,  the  subject  methods  are  all  suggested  for  sub- 
tasks in  defining  the  geoid.  The  most  basic  methods  come  directly  from 
potential  theory:  the  surface  integrals  of  anomalous  potential  (equivalent 
to  undulation)  or  gravity  anomaly.  These  integrals  map  the  anomalous 
potential  or  gravity  anomaly  over  some  closed  surface  onto  the  gravity 
disturbance  vector  at  any  point  outside  the  closed  surface  (see  Section 
H below).  These  integrals  are  approxim'ted  by  partitioning  this  refer- 
ence surface  into  a finite  number  of  elemental  areas  and  forming  the 
sum  of  the  products  of  the  approximated  integrand  and  the  respective 
elemental  area.  The  global  nature  of  this  task  cen  be  ameliorated  by  a 
variable  grid  spacing:  a fine  grid  in  the  immediate  vicinity  of  the  evalu- 
ation nested  in  a sequence  of  grids  of  ever-increasing  coarseness  (Ref 
14:120).  The  number  of  elements  for  a coarse  global  5°X5°  partition  is 
over  2500,  so  the  number  of  parameters  for  such  a model  could  easily 
reach  10,  000.  The  number  of  multiply  and  add  operations  would  be  on 
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that  order  also  for  each  gravity  evaluation. 


Morrison  (Ref  36)  suggests  another  integral  approach  based  on 
integrating  over  a geoidal  surface  density  model.  In  the  geopotential 
application  he  suggests  partitioning  the  geoid  surface  into  1640  equal- 
area  elements.  Although  more  flexibility  can  be  formulated,  the  compu- 
tationally efficient  form  is  to  assume  constant  density  layers  within  each 
block.  In  this  form,  the  model  seems  no  different  for  our  application 
than  the  point  mass  model. 

The  point  mass  model  is  the  prime  candidate  from  the  mass 
distribution  modeling  area.  Other  mass  distribution  models  appear  in 
geophysical  prospecting  (Refs  18  and  37),  but  they  offer  no  special 
advantage  for  our  purpose.  The  mass  distribution  techniques  allow  a 
direct  gravity  calculation  using  Newtonian  gravity  formula,  in  turn,  on 
each  mass  element.  We  can  model  the  anomalous  field  as  closely  as  we 
wish  by  increasing  the  number  of  elements  and  decreasing  the  grid  spac- 
ing (Ref  38).  Such  a modeling  technique  has  global  (Ref  39)  as  well  as 
local  (Refs  24  and  30)  possibilities. 

The  MINUTEMAN  Launch  Region  Gravity  Model  is  an  applica- 
tion of  point  mass  modeling  to  aid  INS  performance.  As  mentioned 
previously,  this  model  was  not  stored  in  and  executed  by  the  airborne 
computer;  the  effect  of  anomalous  gravity  was  compensated  for  in  pre- 
mission targeting  calculations  using  a larger  ground-based  computer. 
The  point  mass  grid  spacing,  similar  to  our  integral  approximation,  is 
based  in  a nested  sequence  of  grids  of  ever-increasing  coarseness. 
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The  finest  grid  is,  naturally,  in  the  area  immediately  surrounding  the 
silo.  The  point  masses  are  submerged  below  the  reference  surface  a 
depth  equal  to  the  grid  spacing  to  enhance  parameter  identification  con- 
vergence (Ref  40:5-6).  This  method  has  built-in  upward  continuation  in 
the  Newton  gravitation  inverse-square  equation.  The  problems  associ- 
ated with  parameter  identification  and  the  computational  burden  of  the 
inverse- square  calculation  for  a large  number  (2520  for  MINUTElvlAN) 
of  point  masses  must  be  considered  when  evaluating  this  modeling  tech- 
nique. A method  which  eliminates  these  costly  global  calculations  might 
prove  more  useful. 

The  local  functional  expansions  might  fill  this  expectation. 
Paraphrasing  Junkins  (Ref  41),  we  may  build  a global  (or  less)  family 
of  locally  valid  functional  expansions  rather  than  one  globally  valid 
aeries  expansion.  These  techniques  are  closely  associated  with  interpo- 
lation techniques --indeed  it  can  be  argued  that  that  is  all  they  really  are. 
Junkins  (Ref  41,  42  and  43)  builds  a general  technique  of  partitioning  the 
shell  of  space  out  to  some  radius  above  the  geoid  into  prisms.  Then, 
gravity  is  modeled  by  a functional  expansion  within  each  prism.  Special 
interpolation  techniques  are  discussed  should  some  order  of  continuity 
be  desired  from  block  to  block.  The  functional  form  is  discussed  in 
general,  but  Chebychev  polynomials  are  stressed.  Other  potential 
modeling  functions  could  come  from  bicubic  spline  functions  (Ref  44), 
multiquadric  equations  (Ref  45),  binary  sampling  functions  (Ref  46),  or 
Walsh  functions  (Ref  46).  Some  extension,  say  from  bicubic  to  tricubic 
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splines,  and  some  additional  development  is  required  to  put  any  of  these 
latter  functicas  in  a form  compatible  with  airborne  use. 

3.  Transformed  Integrals  - 

The  concept  of  applying  integral  transform  techniques  to 
gravity  data  processing  is  not  a new  one.  Geophysical  interpretation  by 
"wave-number"  filtering  techniques  have  been  used  by  the  petroleum 
industry  since  about  1955  (Ref  18:158).  An  example  application  is  solv- 
ing the  inverse  gravity  problem  (mass  distribution  from  gravity  meas- 
urements) to  surmise  the  shape  and  location  of  ore  deposits  (Ref  18: 
179-185).  This  process  requires  a downward  continuation  of  measured 
gravity;  we  are  concerned  with  an  upward  continuation  of  this  same  type 
data. 

Transform  techniques  have  only  recently  been  considered  for 
gravity  survey  purposes.  Heiskanen  and  Moritz  make  no  mention  of 
this  possibility  in  Physical  Geodesy  (Ref  14)  published  in  1967  and  the 
standard  reference  in  the  field.  In  1974,  Davis,  et  al  (Ref  47)  used 
Fourier  transform  analyses  in  comparing  relative  errors  for  several 
algorithms  used  in  computing  vertical  deflections.  Then,  Davis  (Ref  48) 
used  one-  and  two-dimensional  Fourier  transform  error  analyses  as  a 
basis  for  designing  geophysical  surveys.  In  1975,  Long  (Ref  49:44-45) 
suggested  applying  Fast  Fourier  Transform  (FFT)  techniques  to  solu- 
tions of  Stokes  and  Vening-Meinesz  integrals.  More  recently  (1976), 
Thomas  and  Heller  (Ref  50:Chapters  3 and  4)  proposed  a comprehensive 


gravity  data  processing  system  based  on  frequency  domain  techniques. 

These  works  and  suggestions  have  dwelled  on  surface  gravity 
calculations.  The  obvious  extension  is  to  apply  these  methods  to  air- 
borne gravity  calculations.  Prado  (Refs  51  and  52)  has  developed  this 
strategy  using  Hilbert  transforms  to  convert  the  spatial  convolution 
integral  equations  (upward  continuation  and  Vening-Meinesz)  into  spatial 
frequency  domain  multiplications.  Closed-form  expressions  are  pre- 
sented for  the  flat-Earth  case;  the  spherical-Earth  case  remains  an 
area  of  active  research.  Reference  52  also  provides  an  analytical 
approach  for  specifying  the  density  and  extent  of  the  survey  required. 
Presumably,  the  transformed  data  from  a gridded  survey  would  be 
stored  in-flight,  so  gravity  model  parameter  storage  requirements  can 
also  be  assessed. 

A definite  consideration  is  the  existence  of  specific  hardware 
to  p.^rform  the  "butterfly"  operation  (Ref  53:296-297)  that  is  the  heart  of 
the  EFT.  One  can  conceive  of  an  anomalous  gravity  computer  as  a 
separate  functional  module.  Such  a unit  would  accept  navigation  position 
estimates  from  the  airborne  computer  as  input;  would  perform  the 
necessary  transform  inverse  and  interpolation;  and  would  provide 
anomalous  gravity  as  the  output.  Such  a unit  would  allow  this  technique 
to  be  incorporated  in  present  systems  with  modest  interface  and  compu- 
tational burden  on  existing  airborne  computers. 
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4.  Gradiometry 


We  cannot  overlook  the  one  development  which  treats  the 
gravity  error  not  as  a modeling  problem  but  as  a measurement  problem. 
The  gradiometer  research  and  development  addresses  the  real-time 
measurement  of  the  gravity  tensor,  r(r_®),  which  will  allow  computation 
of  anomalous  gravity  in  a manner  which  permits  real-time  compensation 
for  its  effect  on  navigation  estimates.  Although  gradiometers  date  back 
to  the  late  nineteenth  century  experiments  of  Baron  Von  Eotvos,  research 
has  been  concentrated  in  the  last  decade.  The  motivation  for  this 
research  is  the  desire  to  mobilize  the  gravity  survey. 

Gravimeters  used  in  static  gravity  measurements  are  accel- 
erometers which,  in  that  role,  measure  the  acceleration  required  to 
keep  a test  mass  stationary  with  respect  to  an  Earth-bound  observer. 

If  we  use  a gravimeter  in  the  dynamic  environment  of  a mobile  survey, 
say  airborne,  we  must  compensate  for  base  motions.  From  the 
Principle  of  Equivalence  (Ref  1:2),  we  cannot  measure  gravitation 
directly;  so,  the  gravimeter  has  the  same  gravitation  observation  as 
the  accelerometers  in  an  INS.  The  combination  of  an  INS  and  a gradi- 
ometer can  provide  the  basis  for  a statistical  estimate  of  gravitation 
(Refs  54  and  55),  but  these  procedures  are  not  compatible  with  real- 
time compensation.  Moritz  (Ref  56)  demonstrated  that  the  spatial  deriv- 
ative of  gravitation  can  be  measured,  in  principle,  in  a dynamic  environ- 
ment. Conceptually,  we  could  calculate  gravity  from  this  measurement 
through  a spatial  integration  if  we  know  our  path  through  space  and  we 
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are  given  an  initial  condition  for  gravity  at  our  original  coordinates. 

Thus,  a gravity  survey  could  be  conducted.  In  reality,  we  would  need 
real-time  position  information  which  could  come  from  an  INS.  So,  the 
optimal  combination  of  gradiometer  and  INS  evolves  naturally.  Since 
the  INS  needs  gravitation,  either  model  or  measured,  as  am  input,  we 
can  improve  overall  system  performance  by  using  the  results  of  the 
gradiometer  measurements.  To  avoid  the  open-loop  propagation  of 
gradiometer  errors,  the  reference  field  can  be  used  to  make  gradiometer 
biases  observable  {Ref  57  or  Ref  58).  V/'e  need  a mathematical  formula- 
tion for  computing  anomalous  gravity  from  gradiometer  measurements 
and  from  the  reference  field  properties. 

The  first  step  in  this  formulation  is  to  recall  the  definition 

Nbw  consider  the  time  derivative  of  8^  from  an  e-frame  observer's 
point  of  view.  Operating  on  (5)  we  get 


* rm(£®)  - r(r®) 


le 


(23) 


where  the  subscript  e denotes  the  time  derivative  with  respect  to  an 
e-frame  observer.  Extending  (23)  to  other  coordinate  .frames  follows 
from  a straightforward  application  of  vector  calculus.  We  have  ), 


35 


■where  the  subscript  here  refers  to  the  model  not  the  frame,  from  our 
low-order  reference  model.  We  are  proposing  to  measure  r'(*)  with 
gradiometers.  We  can  form  an  estimate  of  dS^/dtg  with  the  additional 
position  and  velocity  from  our  navigation  filter: 

[d's^®)/dt]g  = -r(r®)]  (24) 

Given  initial  conditions,  we  can  estimate  5^  by  integrating  (24). 
/\ 

Given  this  5g,  we  can  estimate  total  gravitation  by  inverting  (5); 


This  demonstrates  the  possibility  of  using  gravity  gradient  measure- 
ments to  calculate  total  gravitation;  again,  we  must  put  such  a calcula- 
tion into  a total  navigation  algorithm  which  will  identify  and  account  for 
gradiometer  bias  for  practical  use  of  the  new  information  (Refs  57  and 
58). 

We  have  discussed  the  use  of  these  gravity  gradient  measures 
without  discussing  how  such  measures  could  be  made.  The  Principle  of 
Equivalence  eliminates  the  accelerometer,  or  gravimeter,  as  a gravi- 
tation sensor.  However  two  accelerometers  in  the  same  dynamic  environ- 
ment, with  input  axes  parallel  and  separated  by  a small  distance  can  be 
used  to  sense  differential  acceleration.  Since  the  dynamics  are  essen- 
tially the  same  for  such  an  accelerometer  pair  mounted  on  a space- 
stable  inertial  platform,  any  differential  acceletation  can  be  attributed 
to  variations  in  gravitation.  We  do  not  require  an  inertial  platform;  if 
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the  accelerometer  pair  is  allowed  to  rotate  with  respect  to  inertial  space 
we  must  sense  this  rotation  and  compensate  for  its  affect  when  process- 
ing the  accelerometer  measurements,  however.  A simple  way  of  seeing 
the  nature  of  the  measurement,  is  to  view  the  two  accelerometers  as  the 
simple  test  masses  they  embody.  One  gradiometer  design  by  The 
Charles  Stark  Draper  Laboratory  is  based  on  this  simple  mass  dipole* 
concept.  If  we  treat  these  test  masses  as  point  masse  > separated  by  a 
lever  arm,  the  difference  in  gravitation  between  the  pa’r  will  create  a 
torque  requirement  based  on  maintaining  a constant  relative  orientation 
of  the  axis  passing  through  the  mass  elements.  This  torque  has  two 
degrees  of  freedom  in  that  a two  dimensional  vector  space  of  torques  is 
required  to  counteract  any  possible  torque  generated  by  gravity  vari- 
ations. Decomposing  this  torque  into  coordinates  perpendicular  to  the 
dipole  axis  and  dividing  each  component  by  the  product  of  elemental  mass 
and  mass  separation  distance  squared  yields  a discrete  approximation  of 
two  components  of  the  gravity  tensor  {see  Figure  2 for  an  example). 

With  this  method  for  computing  the  measured  gravity  tensor 
components,  let  us  investigate  the  number  of  gradiomelers  required  to 
completely  specify  the  full  tensor.  In  Appendix  A we  show  that  the  nine 
elements  of  the  gravity  tensor  are  related  through  LaPlace's  equation 
and  continuity  such  that  only  five  of  the  components  are  independent. 

Then  three  of  our  two-degree-of-freedom  gradiometers  can  provide 

% 

Dipole  is  used  here  for  two  positive  mass  units  which  is  at 
variance  with  the  electromagnetic  use  of  this  term. 
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Fig.  2.  Example  Gradiometer  Output 


measurements  which  span  this  five  dimensional  space--similar  to  two 
two-degree-of  freedom  gyroscopes  spanning  the  three  dimensional 
angular  rotation  space  tor  an  INS. 

With  the  implementation  concepts  outlined,  we  return  to  the 
question  of  gradiometer  errors.  We  have  indicated,  in  general  terms, 
how  the  reference  field  can  be  used  as  an  aid  to  eliminate  gradiometer 
bias  effects  from  our  estimates.  The  gradiometer  error  models  are  in 
a formative  stage.  Some  parametric  studies  have  been  completed  (Refs 
57  through  6l)  which  give  an  indication  of  how  much  relief  we  can  expect 
from  these  instruments.  The  goal  for  these  residual  errors  is  on  the 
order  of  0.  1 EotvSs  units  (1  Ebtvos  unit  = 1 EU  = 10“^sec"^).  This 
goal  is  based  on  gravity  survey,  not  navigation,  requirements.  These 
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studies  indicate  that  substantial  navigation  improvements  are  attained 
■with  instrument  errors  an  order  of  magnitude  higher. 

With  such  devices  in  prototype  testing,  one  might  question  the 
need  for  improving  the  traditional  gravity  modeling  techniques.  The 
fact  is  that  operational  gradiometers  are  no  near-term  certainty.  The 
additional  weight,  space  and  power  required  to  suspend  this  gradiometer 
triad,  in  a manner  which  allows  us  to  isolate  from  or  compensate  for 
inertial  effects  of  rotation,  will  limit  grrdiometer  use  to  relative  large 
systems.  So  gradiometers  are  rot  a panacea  for  our  gravity  modeling 
problems.  They  may,  however,  provide  the  only  economical  method 
for  collecting  the  gravity  survey  data  which  will  support  the  model 
improvements  we  may  propose. 

H.  Basic  Theory 

Each  of  these  candidate  gravity  modeling  improvements  con- 
sists of  a process  for  producing  a local  gravity  estimate  using  available 
gravity  measurements.  The  gravity  data  is  collected  where  and  how 
economical  measurements  can  be  made;  we  require  gravity  estimates 
throughout  the  region  of  possible  INS  operation.  This  practical  problem 
must  be  approached  by  appealing  to  Newtonian  gravitational  theory.  This 
theory  has  developed  over  the  centuries,  but  only  a small  subset  applies 
to  our  problem.  While  some  theoretical  aspects  have  arisen  in  previous 
discussions,  our  attention  was  directed  elsewhere  and  the  theoretical 
aspects  were  incompletely  covered.  The  following  account  recapitulates 
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and  amplifies  those  areas  directly  applicable  to  the  proposed  research. 

In  Newtonian  gravitational  theory,  the  fundamental  information 
is  encoded  in  the  mass  distribution.  We  conceive  an  Earth  mass  distri- 
bution fiinction  expressed  as  a density  (mass  per  unit  volume)  function  of 
the  radius  vector;  p(r).  This  information  is  transmitted  in  the  form  of 
gravitational  force  by  the  relationship 

G(£)  =///  \ (r*  - r)  dw  (26) 

E ir‘  - r 1^ 


where  the  triple  integral  is  over  the  set  E which  encompasses  all  Earth 
mass,  K is  the  Newtonian  gravitational  constant,  £'  is  the  Earth-relative 
radius  vector  to  the  incremental  volume  element  symbolically  called  du, 
and  £ is  the  radius  vector  to  the  point  P where  gravitational  force  evalu- 
ation is  made.  The  gravitational  potential  function,  V(*)»  summarizes 
this  information  as  a scalar  by 


Kp(r») 
|r'  - rl 


du 


(27) 


*We  drop  the  superscript  e since  these  equations  do  not  involve 
time  derivatives  and  are  generally  valid  for  any  choice  of  coordinate 
frame. 
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These  fxmdamental  gr;-.'  itational  quantities  are  related  by 


G^(r)  =dV(r)/dr  (28) 

We  may  derive  (26)  by  formally  applying  (28)  to  (27). 

These  equations  are  necessary  theoretical  concepts,  but  a 
direct  approximation  of  either  (26)  or  (27)  requires  a measure  of  density 
throughout  the  Earth's  volume  which  is  impractical  if  not  impossible. 
Fortunately,  the  information  necessary  to  reproduce  the  external  gravity 
field  is  completely  summarized  in  the  form  of  either  potential  or  normal 
gravitational  force,  Gj^,  over  a closed  surface*  S which  encloses  the  set 
E.  There  are  some  drawbacks  to  this  simplification;  however  we  now 
have  a problem  that  is  merely  impractical,  not  impossible.  That  is,  we 
can  collect  sufficient  Earth  surface  gra.^ity  data  to  produce  useful  repre- 
sentation of  the  Earth's  gravitational  field. 

We  can  measure  gravity  directly  on  the  Earth's  surface  with 
accelerometers,  called  gravimeters,  which  are  held  fixed  with  respect 
to  the  rotating  Earth.  These  measurements  can  be  processed  ("reduced" 
in  the  parlance  of  Physical  Geodesy)  to  produce  a representative  free- 
air  gravity  on  the  geoid  such  that  Earth  mass  above  the  geoid  is  accounted 
for  and  can  be  neglected  (Ref  14;Chapter  3 and  242).  The  gravitational 

^Potential  theory  includes  much  broader  classes  of  closed 
surfaces,  but  this  depth  will  suffice  for  our  purposes.  The  simply- 
closed  surface  is  analogous  to  the  Jordan  curve  from  complex  variable 
theory.  It  is  a two-dimensional,  connected  and  closed  manifold  of  finite 
area  and  enclosing  a finite,  non-zero  volume  (open  set)  in  Euclidian 
three  space. 


potential  can  also  be  measured.  This  indirect  measurement  treats  the 
aea-surfaces  as  (on  the  average)  equipotential  surfaces  and  uses 
aatellite-to-sea  altimetry  as  a measure,  albeit  noisy,  of  the  geoidal 
height  or  undulation,  N,  above  the  reference  surface.  For  example,  if 
the  satellite  altimeter  measure  is  h and  if  we  estimate  the  altitude  using 

A 

both  satellite  ephemeris  and  our  Earth  surface  model  to  be  h,  we  may 
define  a measure  of  geoidal  height  by 

/w  A 

N = h - h.  (30) 


We  get  anomalous  potential  T by  a simple  first  order  calculation 


^ A 

T = 7 N 


(31) 


where  7 is  the  reference  gravity  value  on  the  reference  surface  directly 
below  the  sea  surface  point  representative  of  the  footprint  illuminated 
in  the  measurement.  Gravitational  potential  V varies  from  the  refer- 
ence value  by  this  anomalous  amount  T,  so  we  can  form  a measure  of  V 
by 


A 

V = T + V 


A ^ A ^ 

ref  - (h  - h)  + Vj.g£. 


(32) 


Other  methods  exist  for  measuring  gravity  anomalous  behavior.  Most 
are  based  on  variational  techniques;  for  example,  attributing  satellite 
jrbit  perturbations  to  unmodeled  gravitational  effects  (Ref  14:341-357). 

We  are  concerned  here  with  demonstrating  the  possibility  of 
such  measurements,  not  an  exhaustive  recounting  of  the  methods.  Now, 
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having  demonstrated  the  possibility  of  measuring  either  gravity  or 
potential  over  the  surface  of  the  Earth,  we  can  direct  our  attention  to 
Newtonian  potential  theory.  This  theory  justifies  these  data  sets  as 
sufficient  to  predict  gravitation  throughout  space  above  the  Earth's 
surface. 


We  start  with  the  fact  that  V from  (27)  satisfies  Poisson's 


Equation 


AV  = - 4irKP 


(33) 


where  A is  the  LaPlacian  operator.  In  regions  outside  the  Earth's 
mass  (ignoring  the  atmosphere  and  extraterrestrial  bodies  for  now)  p 
is  zero  so  V satisfies  LaPlace's  Equation 


AV  = 0.  (34) 

Equation  (34)  is  an  elliptic  partial  differential  equation  which  leads  to 
Fredholm  type  integral  equations  (Ref  62:230).  Solutions  to  (34)  are 
called  harmonic  functions  and,  in  general,  form  an  infinite  dimensional 
algebraic  vector  space.  The  general  solution  can  be  expressed  as  a 
linear  combination  of  some  basis  set  of  functions  with  component  coef- 
ficients to  be  determined  by  constraints  placed  on  the  solution  by  the 
particular  problem.  By  selecting  coefficients  which  meet  the  boundary 
conditions  of  our  surface  gravity  measurements,  for  example,  we 

)((  2 a 2 

In  Cartesian  coordinates  A = + 2 — 

dx2  ey2 


produce  a solution  which  both  satisfies  (34)  and  reproduces  the  observed 
data.  The  only  questions  of  such  a solution  are  existence,  uniqueness, 
and  continuous  dependence  on  the  data. 

The  continuity  question  is  easily  answered  since  the  data  is 
inside  the  Fredholm  integral.  As  long  as  the  data  is  reasonably  well 
behaved  (and  we  have  it  completely),  the  derivatives  are  assured.  The 
solution  exists  according  to  DiricHet  Principle  and  is  uniquely  defined 
by  boundary  conditions  of  potential  or  normal  gravity  by  Stoke' s Theorem 
(Ref  14:16-17).  When  the  potential  is  the  boundary  condition,  we  call 
this  the  Dirichlet  Problem;  when  the  normal  gravity  is  the  boundary  con- 
dition we  call  this  the  Neumeinn  Problem. 

Thus,  knowledge  of  V or  Gj^  over  an  enclosing  surface  is  suf- 
ficient information  to  recreate  the  gravitational  field  of  the  Earth's 
mass  distribution.  While  this  condensed  information  can  recreate  the 
gravitational  field,  it  does  not  uniquely  characterize  the  mass  distribu- 
tion which  created  it.  This  non-uniqueness  for  the  inverse  problem- - 
determining  mass  distribution  from  potential  on  a surface- -is  the  bane 
of  geophysical  prospecting  and  Physical  Geodesy  where  the  mass  distri- 
bution is  inherently  valuable  information.  This  issue  is  important  to  us 
since  some  modeling  techniques  (e.  g.  point  mass  models)  are  based  on 
mass  distribution  representation. 

Geophysical,  or  geological,  prospecting  has  been  the  prime 
motivation  for  development  of  the  inverse  solution  techniques.  The 
problem  is  not  only  ambiguous;  it  can  be  ill-posed  as  well  (Ref  40)  since 
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•we  will  never  have  data  at  every  point  on  an  enclosing  surface.  These 
difficulties  have  not  prevented  the  use  of  inverse  techniques  {Refs  18  and 
37)  since  human  interpretation  can  be  used  to  filter  out  obvious  trends. 
Also,  the  recognition  that  most  calculations  czui  be  posed  as  multi- 
dimensional digital  filtering  problems  has  resulted  in  the  more  success- 
ful techniques  (Ref  18:157-185). 

While  we  must  be  aware  of  the  possible  ill-posed  nature  of 
identifying  model  parameters,  we  must  keep  in  mind  that  our  final 
objective  is  to  model  the  gravitational  force,  not  the  underlying  mass 
distribution.  Lee  (Ref  63)  expresses  the  appropriate  point  of  view  for 
this  situation:  objective-oriented  identification.  Interpreting  this  con- 
cept in  terms  of  the  model  identification  tasks  we  shall  face:  We  must 
judge  our  parameter  identification  not  by  how  well  the  Earth's  mass 
distribution  is  represented  but  by  how  well  measured  gravity  data  is 
recreated  and  how  consistent  this  process  is  with  respect  to  gravitational 
theory  extant.  Since  we  are  modeling  for  an  INS,  we  should  judge  the 
final  model  performance  weighted  by  the  spectral  response  of  the  navi- 
gation filter.  That  is,  we  shall  search  for  models  that  are  most  accur- 
ate in  the  passband  of  the  INS  algorithm  witli  special  emphasis  on 
behavior  near  the  Schuler  frequency. 

Returning  now  to  our  discussion  of  (34),  we  may  formally  solve 
for  the  potential  at  any  point  outside  S,  the  enclosing  surface,  by  identi- 
fying the  coefficients  for  the  previously  mentioned  basis  functions.  In 
gravitational  theory,  the  most  common  set  of  basis  functions  are  the 
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spherical  harmonics.  We  generate  these  solutions  by  first  expressing 
(34)  in  spherical  coordinates: 


AV  = + I.  + _1_  _ tan  (f>  d V ^ 1 d^Y  _ q 

d r^  r 0r  0(^)2  0 r^cos^^  dX2 

(34a) 


■where  r is  the  radius,  <{>  is  geocentric  latitude,  and  X is  longitude  (either 
inertial  or  Earth-relative  longitude  applies  in  this  case).  By  separation 
of  variables,  the  solution  Y{x,<i>,\)  is  formed.  The  familiar  spherical 
harmonics  are  the  outcome  of  this  development.  Heiskanen  and  Moritz 
(Ref  14:18-35)  give  a lucid  description  of  this  method.  The  general  solu- 
tion can  be  written  as 


V(r,  <^>,X) 


where  P„  „ is  the  Legendre  function  of  degree  n and  order  m,  and  where 

XI  f XXX 

a and  b_  „ are  the  associated  undetermined  coefficients.  As 
n,  m n,  m 

previously  mentioned  (page  28),  determination  of  these  coefficients  in  a 
global  sense  is  zin  inefficient  means  of  producing  a detailed,  localized 
gravity  model. 

Another  approach  to  solving  (34)  subject  to  boundary  conditions 
on  V or  Gjj  over  S is  to  form  the  associated  integral  equation.  This 
approach  is  simpler  when  the  enclosing  surface  is  a spherical  shell. 
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Our  surface  measurements  can  be  reduced  to  the  geoid  which  is  well- 


approximated  by  a Bjerhammer  sphere.  This  approximation  introduces 
errors  (Ref  14:241-242)  on  the  order  of  Earth  flattening  (1/298.  26,  Ref 
17:16)  which  is  troublesome  if  we  dead  with  full-valued  V and  This 

problem  is  ameliorated  if  we  subtract  the  reference  field  contribution 
from  the  reduced  data.  The  remaining  anomalous  quantities  can  be 
treated  as  first-order  errors  and  since  the  flattening  is  of  the  same 
order,  the  resulting  products  can  be  neglected  as  second  order  terms. 
This  method  is  analogous  to  INS  modeling  techniques  (Ref  1),  and  thus, 
this  approximation  of  the  geoid  by  a sphere  is  consistent  with  our  other 
analytical  tools . 

This  approach  leads  to  the  Poisson  Integral  characteristic  of 
the  upward  continuation  of  a harmonic  function: 

T(r,4>,  X)  = f fz  A')  cos  4>'d(j>'dX'  (36) 

4ir  J J 

X=:0  <{)*=-£ 
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where  R is  Bjerhammar  sphere  radius;  is  the  radius  vector  to  the 

sphere  surface  incremental  area  defined  by  geocentric  latitude  4>'  and 

longitude  \';  and  r_  is  radius  vector  to  point  defined  by  r, 4>»  and  X (Ref 

14:37).  Hereinafter  the  set  |r,4>,\}  and£will  be  used  interchangeably; 

also,  the  integral  over  the  unit  sphere  will  be  abbreviated  to/  ( ) dcr 

•'cr 

where  <r  represents  unit  sphere  surface  (the  limits  of  integration  in 
(36))  and  dcr  represents  cos  4>'d«()'dX’. 


This  equation  assumes  we  have  as  input  data  the  function  T(^') 
over  the  surface  of  the  reference  sphere.  Another  form  for  expressing 
T(^')  in  integral  form  comes  from  the  relationship  of  T to  the  gravity 
anomaly  Ag  (Ref  14:89): 

Ag  = - 8T/dr  - 2T/r  . (37) 

This  relationship,  applied  to  (36),  leads  to  Stoke 's  Formula: 

T(r)  = Ag(R')  S(r,*)  da  (38) 


■where  is  the  central  angle  between  r^  and  and  where  S(r,  >1^)  is  Stoke's 
function  given  by  (Ref  14:233) 


S(r,^)  = _M_  + R ,3RIR«-rl  . M cos  ^ 


IR'-rl  r 


^2  ^2 

r r 


[5  + 3 In  ^r-R  t IR-rij 


(39) 


Equation  (38)  requires  input  data  of  gravity  anomaly  over  the 
sphere.  The  set  (36)  and  (38)  then  provide  two  direct  methods  of  com- 
puting T(r)  from  measured  data.  Since  we  ultimately  want  the  gravity 
disturbance  vector  function,  5£(r),  we  can  use 


fi£(r)  = dT/dr  . 


(40) 


For  spherical  coordinates,  this  becomes 
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«gj=  dT/ar, 


(40a) 


and 


6g.  =i.  3T/a«  , 
« r 


(40b) 


1 

r cos  ^ 


dT/ax. 


(40c) 


We  can  now  apply  (40)  to  (36)  and  (38),  in  turn,  to  yield  integral  relation- 
ships between  the  input  data  (T  or  Ag,  respectively)  and  the  desired 
final  result  In  applying  (40)  it  is  convenient  to  express  ^ in  terms 
of  the  latitudes  and  longitudes  or  the  defining  vectors  and  r^: 


'll  = cos"^  [sin  <t>  sin  0'  + cos  0 cos  0'  cos  (X'  - X )]  , (41) 


Now  since  iR'-rl  = !R^  + r^  - ZrRcos 0 1 


(42) 


by  the  Law  of  Cosines,  v/e  have 

17?  , .1/2 
= I R'^  + r - ZrR  [sin  0 sin  0'  + cos  0 cos  0 ’ cos  (X*  - X)] ) 

(43) 

With  (43)  the  necessary  partial  derivatives  are  straightforward,  and  we 
get  from  (36); 


5g^  = ^ /* M(r,0)T(R')dc 

^ 4» 


(44a) 


where 


^ 5R^r  - r^  - Rr^cos  0 - 3r3  cos  0) 

M(r,  0)  = z ,(Ref  14;37) 

JR>  - r I ^ 


(44b) 


_ 3R^(r^-R^)  /*  T(^')  [cos  0 sin0'  - sin^cos  0'  cos  (>'  -X)]  do 

iiTTp 

(44c) 

I 

e_  _ 3R2(r2-R2)  /*  t(R»)  cos  0‘  sin  (X*  - \) 

'«X j; 7 -=^— 5 do  • (44d) 

o iR'  - r I ^ 


Now,  in  applying  (40)  to  (38)  it  is  convenient  to  define  an  azimuth  a of 
the  line  segment,  or  arc,  fromir  to^'  (see  Figure  3). 


tan  « = S.9£-»*  (X*  - X) 

cos  0 sin0'  - sin  0COS  0 'cos  (X ' - X) 

And  note  that  80/d0  = - cos  a 


and 


30 /dX  = - cos  0 nin  ft  . 


Then,  (38)  becomes 


dS(r,0) 

3r 


da. 


(45) 

(46a) 

(46b) 


(47a) 


cos  a da, 

'*x  = - f;  cino  d. 


where 


(47b) 


(47c) 


(47d) 


presented  in  Reference  14  on  page  234. 

Equations  (44)  and  (47)  present  the  direct  methods  for  computing 
the  disturbance  vector  from  either  type  of  input  data:  T or  Ag*  Since 
this  disturbance  vector  represents  the  error  in  our  model  or  reference 
gravity  formula,  we  may  compensate  the  reference  model  by  adding  this 
computed  disturbance  acceleration  to  the  reference  gravitation  accelera- 
tion vector.  Hence,  our  algorithm  for  approximating  either  (44)  or  (47^ 
provides  one  alternative  gravitational  modeling  improvement  as  we 
discussed  in  the  last  section. 

Another  method  for  calculating  the  disturbance  vector  uses  the 
technique  of  replacing  the  disturbing  rriasses  by  an  equivalent  surface 
mass  layer  (see  Theorem  of  Chasles  (Ref  14:13). Once  this  surface  mass 
distribution  function  has  been  defined  from  the  input  data,  we  simply 
apply  (26)  taking  advantage  of  the  fact  that  the  volume  integral  can  now 
be  replaced  by  a surface  one.  This  "coating  method"  requires  both  types 
of  input  data  (T  or  N and  Ag)  to  compute  the  surface  density  function 
(Ref  14:236-238).  The  resulting  formulae  are  simpler  than  those  for  the 
direct  methods.  This  factor  might  offset  the  costs  of  additional  meas- 
urements. and  calculation  in  forming  the  surface  density  function. 

These  three  methods  (Poisson  Integral,  Stoke's  Integral  and 
coating)  have  approximately  the  same  accuracy  (Ref  14:243).  To  imple- 
ment one  of  these  integral  methods,  we  must  face  the  issue  of  a discrete 
approximation  method.  One  can  see  for  the  coating  method  that  the  most 
irnportcuit  data  is  that  near  ^ = 0.  This  idea  carries  over  to  the  other 
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direct  methods  but  with  less  intuitive  appeal  than  the  surface  distribu- 
tion case.  Hirvonen  and  Moritz  (Ref  64)  considered  this  factor  in  a 
comparison  study  on  the  direct  Stoke's  (47)  and  the  coating  methods. 

The  idea  is  to  investigate  the  effect  of  ignoring  distant  data  in  approxi- 
mating these  integrals.  They  statistically  predict  the  root-mean- 
square  (rms)  residual  gravity  for  each  component  of  the  disturbance 
vector  when  the  integration  set  a is  approximated  by  the  points  within  a 
central  angle  'f'^  of  £.  The  statistical  error  estimation  made  use  of  the 
anomaly  covariance  function  (see  Appendix  C)  as  a basis  for  statistically 
characterizing  the  effects  of  the  distant  zones. 

To  interpret  the  results  of  this  study,  we  should  point  out  that 
each  component  of  the  disturbance  vector  on  the  Earth's  surface  is  an 
approximately  35  mgal  rms  process  (Ref  50).  This  figure  corresponds 
to  a = 0 when  the  compensating  integral  covers  no  area.  For  a 
= 9°»  the  estimated  rms  residual  is  down  to  8 mgal  and  6 mgal  for 
the  direct  and  coating  methods,  respectively.  So,  covering  0.67o  of  the 
Earth's  surface  decreases  our  expected  rms  gravity  disturbance  to,  on 
the  order  of,  20%  of  the  uncompensated  level.  We  require  a central 
angle  coverage  of  90°,  or  one-half  the  Earth's  surface,  to  lower  the 
estimated  rms  disturbance  to  10%  of  the  uncompensated  value  demon- 
strating a severe  Law  of  Diminishing  Return  for  increasing  i^q. 

We  conclude  that  we  will  require  a gravity  survey  over  a vast 
area  to  support  eUiy  of  these  integral  approaches'.  We  can  diminish  the 
computational  burden  by  variable  grid  spacing  of  the  area  necessary  for 
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adequate  compensation  accuracy.  This  variable  spacing  policy  was 
mentioned  earlier  and  is  the  basis  of  point  mass  modeling  techniques. 
We  can  use  the  integral  kernels  and  our  knowledge  of  expected  flight 
path  to  select  a grid  schedule  of  approximately  equal  statistical  influ- 
ence. That  is,  the  error  introduced  by  the  finest  grid  spacing  terms 
should  be  on  the  same  order  as  that  expected  from  the  coarsest  grid 
spaces.  Such  a technique,  it  must  be  noted,  will  place  flight  path 
restrictions  on  the  mission  if  we  intend  tc  salvage  the  accuracy  for 
which  we  are  building  the  compensation. 

The  problem  of  approximating  these  integrals,  as  noted  in  the 
previous  section,  might  be  handled  with  integral  transform  methods. 
The  computational  structure  of  the  Fast-Fourier  Transform  makes  it 
attractive.  The  data  could  conceivably  be  transformed  pre-mission 
leaving  only  the  interpolative  inverse  transform  as  real-time  computa- 
tional burden.  Some  research  is  required  to  adapt  these  FFT  methods 
to  the  variable  grid  spacing  that  we  must  consider. 

These,  then,  are  some  of  the  basic  theoretical  concepts  and 
concerns  relating  to  our  gravity  modeling  task.  We  are  at  a point 
where  concrete  plans  can  be  made  for  future  research.  The  subject  of 
gravity  modeling,  as  we  have  seen,  is  broad  and  reaches  into  the  fields 
of  Physical  Geodesy  and  Inertial  Navigation.  Some  thought  must  be 
given  to  a proper  division  of  responsibilities  in  the  oveiail  resea.. ob  of 
this  subject. 
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■J.  The  Geodesy  Connection 


To  get  the  total  view  of  this  process,  it  is  instructive  to  consider 
a signal  flow  graph  of  the  "information"  stemming  frcm  the  mass  dis- 
tribution and  terminating  in  the  navigation  estimates.  With  some 
poetic  license,  Figure  4 presents  such  a view.  This  figure  graphically 
displays  the  chronological  processing  of  gravitational  information  with 
emphasis  on  the  new  tasks  associated  with  the  augmentation  of  the  refer- 
ence field  model.  We  shall  discuss  the  nature  of  these  tasks  which 
compensate  the  reference  field.  We  shall,  also,  describe  how  the 
design  error  budget  provides  criteria  for  model  performance  and,  in 
general  terms,  how  the  model  selection  process  should  logically  be 
conducted.  These  concepts  are  important  because  the  proposed 
research  will  be  conducted  along  lines  compatible  with  model  selection 
even  though  it  will  not  include  a complete  trade-off  study  for  every  con- 
ceivable model  compensation  method.  Finally,  as  mentioned  previously, 
we  shall  discuss  the  logical  assignment  of  these  new  tasks  to  the  fields 
of  Physical  Geodesy  and  Inertial  Navigation. 

Before  we  cover  these  interface  concerns,  we  need  to  reestablish 
what  the  term  "model"  means.  As  shown  in  Figure  4,  the  complete 
gravitational  model  consists  of  two  parts;  (1)  The  reference  gravita- 
tional field  based  on  a surface -approximating  ellipsoid  and  (2)  the  com- 
pensation model  which  calculates  local  variations.  The  standard  names 
and  notation  become  troublesome  at  this  time  since  the  disturbance 
vector,  which  has  been  treated  as  a perturbation  quantity,  is  now 
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being  approximated.  We  need  an  unambiguous  name  and  symbol  for  the 
new  residual  gravitational  modeling  error.  We  shall  continue  to  call 
the  variation  from  the  reference  field  5^.  The  ellipsoidal  field  in 
earlier  discussions  was  synonymous  with  the  model.  Since  our  new 
model  will  include  this  ellipsoidal  one  as  an  element,  we  shall  refer  to 
it  as  the  reference  field,  symbolized  by  )•  We  shall  refer  to  the 

compensation  model  as  such  and  use  the  symbol  ).  The  remain- 

ing  residual  disturbance  will  be  symbolically  £(• ).  With  these  new 
definitions  (5)  on  page  9 should  be  rewritten  as 

- G ( r ®) . ( 5a) 

Our  complete  model  is  now 

=Gj.gf(r®)  - (48) 

The  traditional  modeling  methods  can  be  viewed  as  a degenerate  form 
of  (48)  with  defined  as  a constant  zero  vector  function.  The 
residual  disturbance  vector  is  then  defined  by 

- G(re)  =-[5£jj,(r.e)  - 6£(re)].  (49) 

In  l''igure  4,  prefers  to  8£^(£)  and  prefers  to^^(£).  With  the 
notcition  issue  resolved,  we  can  discuss  how  this  new  modeling  process 
will  be  accomplished. 

Again  returning  to  Figure  4,  we  will  require  an  extensive  local 
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gravity  survey  to  provide  the  short  spatial  wavelength  information 
necessary  for  model  improvement.  This  survey  data,  taken  over  a 
grid  convenient  for  measurement,  must  be  reduced  to  the  reference 
surface  and  interpolated  to  form  a data  s.et  on  the  grid  which  we  shall 
establish  for  model  parameter  identification.  This  interpolation  will 
undoubtedly  be  necessary  since  the  parameter  identification  algorithm 
must  work  for  arbitrary  survey  grid  structure.  The  parameter  identi- 
fication and  functional  evaluation  blocks  are  shown  for  four  example 
compensation  methods.  We  can  see  the  variety  of  forms  these  tasks 
might  take.  The  parameter  identification  varies  considerably  from 
possibly  using  the  interpolated  grldded  data  directly  as  a parameter  set 
to  the  search  procedures  already  developed  for  point  mass  modeling 
(Ref  30).  After  the  modeling  method  has  been  selected  and  the  param- 
eters have  been  identified,  the  model  is  completely  defined.  The  next 
step,  functional  evaluation  must  occur  in  the  navigation  computer  for 
real  time  compensation  since  the  position  estimate  is  required.  Once 
this  gravitational  disturbance  estimate  is  combined  with  the  estimated 
reference  field,  the  INS  algorithm  is  unchanged.  Figure  4 provides 
examples  of  how  such  compensation  might  be  accomplished.  Clever 
computational  schemes  might  eliminate  some  of  the  blocks  but  the 
equivalent  of  the  four  tasks  must  be  accomplished:  survey,  process  to 
computational  grid,  identify  parameters  and  functionally  evaluate. 

The  new  tasks  fall  into  categories  based  on  the  discipline  primarily 
responsible  for  the  theoretical  development  and  implementation  of  the 


task.  The  gravity  survey  and  subsequent  data  processing  are  in  the 
province  of  Physical  Geodesy.  The  real-time  functional  evaluation  is 
obviously  in  the  field  of  Inertial  Navigation.  The  remaining  tasks  of 
model  selection  and  parameter  identification  are  at  the  intersection  of 
these  fields  of  interest.  The  functional  evaluation  accuracy  is  so 
dependent  on  parameter  identification  criteria  that  these  tasks  right- 
fully belong  together.  The  model  selection  has  a pervasive  effect  on 
the  nature  of  tasks  performed  in  either  area.  The  choice  of  compensa- 
tion method  will  rest  on  a trade-off  which  is  difficult  to  describe  in 
generality.  A specific  design  setting  provides  the  constraints  and  goals 
to  structure  such  a trade-off  without  lapsing  into  abstruse  generalities. 

Let  us  construct  such  a design  scenario  to  understand  how  the  com- 
pensation selection  process  might  logically  be  conducted.  The  navigation 
subsystem  will  typically  be  given  an  accuracy  criteria  based  on  some 
overall  system  performance  requirement.  This  crit  ria  is  based  on  an 
expected  environment  and  is  usually  statistical  in  nature  (e.  g.  rms  posi- 
tion, velocity,  and  attitude  levels  throughout  the  mission).  Suppose  we 
have  been  handed  such  an  error  budget.  Furthermore,  suppose  we  have 
performed  some  analyses  along  the  lines  suggested  by  Levine  and  Gelb 
(Ref  21)  which  indicates  that  the  ellipsoidal  reference  field  is  unaccept- 
able as  a total  model.  Such  a scenario  might  be  a strategic  bomber 
mission  which  precludes  electro-magnetic  emissions  for  detectability 
and  which  assumes  other  navigation  aids  will  ndt  be  available.  We  con- 
clude then  that  we  must  provide  some  compensation  to  the  reference 


field  to  decrease  the  effects  of  gravity  disturbances.  The  INS  error 
budget  must  be  decomposed  to  budgets  for  the  elements  of  the  INS  (e.g. 
accelerometers,  gyroscopes,  gravitation  model,  and  algorithm  numer- 
ical quanti2sation) . With  this  parceling  out  of  the  error  budget,  the 
gravitational  model  will  have  a specific  accuracy  goal  to  meet.  The 
next  step  is  to  select  a compensation  scheme  which  meets  this  derived 
criteria  and  which  is  least  burdensome  in  some  overall  sense. 

"Least  burdensome"  is  subjective  since  many  diverse  costs  must 
be  considered.  The  spectral  content  of  the  gravity  disturbance,  if 
known  or  approximated,  will  allow  us  to  specify  a minimum  spatial 
frequency  to  accomplish  the  accuracy  goal.  By  the  Shannon  sampling 
theorem  this  implies  a minimum  density  for  the  gravity  survey.  The 
extent  of  the  survey  will  depend  on  the  convergence  properties  of  the 
compensation  method.  As  we  have  seen,  the  coating  method  tends  to 
converge  faster  (in  a spherical  cap  size  sense)  than  the  direct  Stoke's 
integral  method.  The  coating  method  requires  two  types  of  input  data, 
however,  so  survey  extent  is  not  a full  account  of  survey  cost. 

The  data  processing  task,  also,  depends  heavily  on  the  compen- 
sation model  choice.  The  integral  techniques  and  the  point  mass  model 
need  data  over  a two-dimensional  surface.  The  functional  approxima- 
tion methods,  however,  need  a three  dimensional  array  of  information 
for  the  three-dimensional  function  parameters  (e.g.  Ref  43).  This 
processing  of  raw  survey  data  and  model  parameter  identification  will 
likely  occur  pre-mission  with  relatively  large  computational  facilities 
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used.  So,  the  costs  may  be  less  onerous  in  this  area  if  the  method  is 
particularly  efficient  in  the  functional  evaluation  phase. 

On  the  implementation  of  a real-time  computation,  the  constraints 
and  costs  are  primarily  associated  with  the  airborne  computer  memory 
storage  and  computational  cycle  time.  The  fine  gravity  model  for 
Reference  24,  for  example,  requires  over  10,000  point  masses.  Each 
point  mass  requires  four  parameters;  so,  over  40,000  memory  loca- 
tions might  be  required  for  parameter  storage  alone.  Furthermore, 
the  functional  evaluation  requires  the  summation  of  the  Newtonian  gravi- 
tational acceleration  from  each  point  mass.  The  impact  of  such  mas- 
sive computational  and  storage  requirements  can  be  put  in  perspective 
when  you  consider  that  operational  flight  computer  memories,  today, 
have  less  than  40,  000  words  for  data  and  program  storage.  Prototype 
computing  systems  with  mega-word  capability  are  available.  So,  future 
systems  may  have  the  flexibility  to  consider  these  more  complete  com- 
pensation models. 

The  crux  of  all  these  costs  is  the  compensation  method.  It  is 
clear  from  the  costs  associated  with  the  survey  and  the  real-time  eval- 
uation that  some  comprehensive,  coordinated  model  selection  study 
must  be  performed  for  each  different  mission  and  system.  While  the 
proposed  research  does  not  cover  the  total  model  selection  issue,  it  is 
important  to  realize  these  issues  to  keep  the  proposed  work  in  context. 
The  structure  of  the  study  and  the  motivation  for  many  assumptions  and 
procedures  will  be  founded  on  these  considerations.  That  is,  we  shall 
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aasume  that  the  least  costly  method  to  reach  a prescribed  mission 
I navigation  accuracy  is  the  driving  force  behind  what  we  do. 

The  proposed  research  will  concentrate  on  the  Inertial  Navigation 
aspects  of  the  problem.  We  shall  assume  that  the  Physical  Geodesy 
I connection  exists  and  will  provide  the  required  gravitational  data  on  a 

computationally  convenient  grid.  Thus,  we  shall  concentrate  on  mini- 
mizing the  real-time  INS  costs.  Again,  the  proposed  research  should 
be  viewed  in  the  context  of  providing  information  to  answer  the  larger 
issue  of  model  compensation  selection. 
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m.  Research  Topic 


The  background  of  INS  gravitational  modeling  shows  the  grow- 
ing need  to  develop  and  demonstrate  techniques  to  compensate  the 
traditional  (usually  ellipsoidal)  INS  gravitational  model.  The  research 
proposed  in  the  following  sections  addresses  the  system  design  prob- 
lem of  selecting  a compensation  modeling  concept  and  of  selecting  the 
level  of  model  complexity  (i.  e.  number  of  j>arameters)  consistent  with 
design  constraints.  The  problem  is  placed  in  this  context  by  address- 
ing the  design  trade-off  between  the  system  performance  cost  and  the 
system  resource  cost  associated  with  the  gravitational  model.  The 
specific  objective  is  to  provide  the  analytical  means  to  evaluate  model 
performance  at  the  system  level  so  that  one  can  choose  both  the  model 
concept  and  the  degree  of  model  complexity.  We  shall  discuss  below 
the  general  setting  for  the  overall  analysis  required  to  select  a gravi- 
tation model.  The  problem  is  partitioned  in  such  a manner  that  the 
fundamental  analysis  required  is  to  determine  INS  performance  for  a 
finite  set  of  design  missions.  All  higher  level  problems  can  be 
approached  given  this  basic  analysis  capability.  An  example  is  sug- 
gested to  demonstrate  the  analysis  method  since  many  of  the  analysis 
steps  depend  on  the  modeling  concept. 


A,  Study  Context 

An  abstract  study  of  gravitational  modeling  errors  seems  a 
barren  exercise  to  one  aware  of  the  strictures  of  an  operational,  real- 
time environment.  The  choice  of  a gravitational  model  is  made  in  a 
system  design  context.  The  gravitational  model  is  one  component  of 
the  INS  which,  in  turn,  is  an  element  of  some  larger  system  designed 
to  perform  some  mission.  The  gravitational  model  is  not  intrinsically 
interesting  from  this  total  system  perspective.  We  are  interested  in 
the  effect  modeling  errors  have  on  system  performance  and  the 
impact  of  the  model  on  such  system  resources  as  computer  memory. 

A conflict  exists  between  our  desire  for  increased  performance  and 
for  efficient  use  of  resources.  An  increase  in  model  complexity  can 
yield  improved  system  performance,  but  it  implies  additional  com- 
puter usage.  This  conflict  generates  a system  design  trade-off  which 
requires  evaluation  of  both  model  performance  costs  and  model 
resource  costs  to  answer  questions  such  as: 

1.  Given  X amount  of  computer  memory  (or  computation  time) 
which  model  concept  at  what  degree  of  complexity  yields  the  best 
system  performance?  and 

2.  Given  Y system  performance  criteria  (say  rms  position 
accuracy)  which  model  requires  the  least  computer  resources? 

The  context  of  this  study  will  be  to  provide  the  analytical  methods  to 
answer  such  design  trade-off  questions.  This  approach  should  make 
the  resulting  methods  immediately  applicable  and  the  motivation 
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more  comprehensible. 


We  shall  concentrate  attention  on  the  performance  cost  evalu- 
ation, but  we  will  need  a measure  of  system  resource  cost  as  a foil 
in  the  trade-off  exercise.  The  modeling  concepts  we  shall  consider 
have  the  quality  that,  barring  other  errors,  they  can  represent  true 
gravitation  to  any  desired  level  of  residual  by  increasing  the  number 
of  model  parameters  (e.g.  more  point  masses).  Model  complexity 
shall  be  used  synonymously  with  number  of  parameters.  The  compu- 
tational overhead  for  each  model  concept  is  the  algorithm  storage 
requirement.  For  a given  model  concept  this  requirement  should  not 
be  strongly  influenced  by  the  number  of  parameters  to  be  processed. 
We  shall  neglect  this  algorithm  storage  cost  in  our  analysis  since  we 
shall  concentrate  on  one  modeling  concept  as  an  example.  Of  course, 
this  cost  must  be  considered  in  any  comparison  between  different 
modeling  concepts  (e.g.  point  mass  versus  Stoke's  Integral). 

The  other  side  of  this  system  design  trade-off  is  the  system 
performance  cost.  We  consider  those  applications  where  a system 
"miss"  is  directly  associated  with  navigation  accuracy.  The  system 
miss  is  assumed  to  be  in  the  form  of  a miss  vector  which  has  com- 
ponents we  wish  to  drive  to  i.ero  (e.g.  downrange  and  crossrange 
miss  for  a ballistic  weapon  delivery).  We  assume  a system  cost 
function  is  associated  with  this  miss  vector.  This  cost  function  maps 
the  miss  vector  space  into  the  real  numbers  and  provides  a means  of 
ordering  miss  vectors  as  either  more  or  less  costly  in  comparison 
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with  other  system  misses.  Our  system  design  must  account  for  a 
wide  range  of  operational  environments  for  a given  mission.  So,  we 
want  to  decrease  the  expected  value  of  this  system  performance  cost 
over  the  range  of  missions  and  environments.  The  nature  of  this  cost 
function  and  the  nature  of  the  statistical  expectation  are  important  in 
formulating  our  analysis  method. 

From  a system  perspective,  we  desire  a measure  of  the  system 
performance  cost  which  reflects  the  evaluation  of  system  cost  function 
for  .every  conceivable  mission  and  environment- -a  weighted  sum  or 
integral.  The  weighting  should  include  not  only  the  probability  of  a 
given  mission  but  the  relative  importance  of  that  mission  with  respect 
to  other  missions,  as  well.  Such  abstract  analysis  is  not  performed 
in  system  design  studies  because  the  resources  to  conduct  such  studies 
are  J,imi.ted.  The  integral  over  all  missions  is  approximated  by  a 
discrete  summation.  The  mission  space  is  subdivided  into  mission 
regions  each  of  which  is  represented  by  one  design  mission. 

The  ICBM  provides  an  easily  understood  example  of  this  mis- 
sion partition  concept.  For  a given  target  range  and  azimuth  the 
ICBM  usually  has  two  solutions  to  the  two  point  boundary  value  prob- 
lem; burn-out  state  vector  to  hit  a specific  target.  This  generates  a 
discrete  (binary  in  this  case)  partition  of  the  mission  space  into  low 
and  high  trajectories.  The  range  and  azimuth  offer  a two-fold  con- 
tinuum for  mission  partitions.  A design  mission  might,  for  example, 
be  specified  as  the  high  trajectory  for  range  X and  azimuth  Y.  This 
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design  mission  would  represent  in  our  system  studies  all  high 
trajectory  missions  over  a set  of  ranges  within  X + AX  and  over  a 
set  of  azimuths  within  Y + AY.  A collection  of  such  design  missions 
can  be  analyzed  and  the  results  combined  with  appropriate  weighting 
to  approximate  the  integral  over  the  mission  space. 

The  mission  space  for  a strategic  bomber,  for  example,  is 
much  more  complex  due  to  the  system's  flexibility.  This  diversity 
would  be  reflected  in  a higher  number  of  design  missions --not  a 
change  in  the  analysis  method.  So,  we  can  focus  our  attention  on  pro- 
viding the  analysis  methods  to  characterize  the  gravitational  model 
system  performance  on  a restrictive  mission  type.  A mission  taxon- 
omy can  be  based  on  geometry  and  the  geography.  The  geometrical 
qualities  one  might  use  are,  for  example,  range,  altitude  and  azimuth 
relative  to  downrange.  If  mission  velocity  is  not  prescribed  by  these 
geometrical  considerations,  the  mission  partition  must  include  these 
variations.  The  mission  type  is,  also,  characterized  by  such  geo- 
graphic considerations  as  the  area  of  mission  origins  and  the  azimuth 
of  the  downrange  grou-'idtrack.  Both  geometric  and  geographic  terms 
cdfect  the  system  level  performance  of  the  gravitational  model. 

The  design  mission  can  be  used  to  characterize  the  effects  of 
anomalous  gravitation  for  a region  of  mission  space.  The  variations 
of  gravitation  within  the  geographic  bounds  of  the  mission  region  is 
the  underlying  random  process.  The  INS  error  dynamics  are  also 
influenced  by  the  geographic  region  (e.g.  Foucault  mode,  see 
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Ref  1:128).  These  error  dynamics  are  also  influenced  by  the  mission 
geometry- -primarily  the  radius  magnitude.  The  mission  velocity 
trcuxslales  the  gravitational  disturbance  quantities  from  the  spatial  to 
temporal  domain.  In  this  manner  the  design  mission  can  be  used  as 
an  analysis  tool. 

To  realize  this  goal  we  must  account  for  the  variation  of  gravi- 
tation within  the  geographic  region  of  interest.  The  gravitational 
disturbance  information  for  this  region  can  be  statistically  summar- 
ized in  terms  of  reference  surface  covariance  functions  of  anomalous 
gravitation  quantities.  We  want  to  include  the  effect  of  gravitation 
over  this  entire  region  since  the  gravitational  disturbance  along  a 
design  mission  trajectory  may  not  be  truly  representative  of  the  entire 
region.  So,  we  shall  seek  the  performance  cost  estimate  for  a given 
region  as  a statistical  estimate  over  all  possible  gravitational  disturb- 
ances in  the  region  and  structure  the  analysis  for  those  reference 
surface  covariance  functions  we  expect  to  be  available  from  survey 
data  (see  Appendix  C). 

The  variations  of  the  INS  error  propagation  model  with  geography 
is  predictable;  therefore,  we  can  partition  the  mission  space  in  such  a 
manner  that  the  design  mission  INS  error  propagation  model  repre- 
sents well  the  entire  mission  region.  A similar  partition  argument 
can  be  made  for  the  geometrical  variations  of  the  mission  space. 
Therefore,  we  shall  assume  the  mission  space  is  partitioned  such  that 
the  INS  error  propagation  along  the  design  mission  is  acceptable  over 
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the  entire  range  of  missions  within  the  mission  region.  This  assump- 
tion is  essential  to  define  a manageable  analysis  problem.  The  vari- 
ations of  INS  error  dynamics  with  geography  are  smooth  enough  that 
the  mission  space  partitions  on  this  account  might  be  widely  spaced. 

The  system  cost  function  must,  also,  be  understood  befox-e 
practical  analyses  can  be  planned.  The  gravitational  modeling  errors 
excite  the  INS  error  state  vectors  through  well-known  (Ref  1)  differ- 
ential equation  models.  The  specific  disturbance  time  function 
depends  on  the  particular  trajectory  (position- time  history).  The  INS 
errors  which  result  from  the  disturbance  input  will  cause  the  overall 
system  to  miss  some  system  target.  A miss  vector  space  is  defined 
in  terms  of  these  component  deviations  from  the  coordinated  system 
objectives.  While  spatial  objectives  are  familiar  (e.g.  downrange 
and  crossrange),  we  can  have  non-spatial  objectives  as  wsll.  An 
error  in  time-on-target  or  attitude  errors  during  antenna  pointing 
operation  are  example  non-spatial  system  miss  quantities.  For  our 
analysis,  we  shall  consider  the  mapping  from  INS  error  state  vector 
space  to  system  miss  vector  space  to  be  linear  by  use  of  well  known 
linear  perturbation  methods.  That  is, 

i(t)  = H(t)X(t)  (50) 

where  ^ is  the  system  miss  vector  corresponding  to  an  INS  error 
state  vector  The  techniques  developed  can  be  generalized  to  a 
broader  class  of  mappings  (i.  e.  convex)  but  the  linear  model  should 
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apply  to  most  problems  of  interest. 

We  assume  the  miss  vector  components  are  each  important  to 
system  effectiveness  such  that  an  increase  in  the  magnitude  of  any 
component  implies,  by  itself,  an  increase  in  the  system  cost  function. 
The  system  cost  fvinction  is  assumed  to  be  a positive  definite,  convex 
functional  over  the  system  miss  vector  space.  This  assumption 
should  certainly  characterize  the  performance  for  a region  near  the 
origin  of  the  miss  vector  space.  The  motivation  for  such  assumptions 
is  to  describe  those  problems  for  which  the  analysis  can  concentrate 
on  the  INS  error  state  statistics.  To  reach  this  goal,  we  must  have 
the  capability  of  reflecting  a system  cost  criteria  back  onto  the  miss 
vector  space  and  then  onto  the  INS  error  state  space.  The  effect  of 
large  scale  miss  on  system  effectiveness  is  not  necessarily  well- 
modeled  by  such  convex  functionals.  The  probability  of  kill  envelopes 
for  nuclear  weapons,  for  example,  define  markedly  non-convex  sets 
due  to  the  various  effects  of  the  weapon.  In  such  cases  a convex  cost 
functional  can  be  derived  by  forming  a performance  index  from  the 
maximum  or  minimum  level  within  a radial  distance  from  the  aim 
point.  The  alternative  is  to  generate  the  statistical  performance  index 
by  Monte  Carlo  analysis  techniques. 

With  the  convex  functional  assumptions,  we  can  formalize  the 
relationship  between  system  cost  criteria  and  the  gravitational  model- 
ing errors.  This  relationship  is  normally  established  through  a 
system  error  budget. 
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The  system  error  budget  is  a design  tool  used  to  distribute  the 
performance  cost  criteria,  through  analysis  and  assumption,  to  vari- 
ous system  error  sources.  Once  established  the  error  budget  forms 
the  subsystem  design  goals  for  a design  iteration.  The  cost  is 
assumed  to  be  stated  as  a level  on  the  expected  cost  (e.g.  circular- 
error-probable  less  than  1 n.mi).  This  level  criteria  yields  an  inter- 
val on  the  positive  half  of  the  real  number  line  which  corresponds  to 
acceptable  expected  system  performance  cost.  The  inverse  image  of 
this  acceptable  interval  onto  the  miss  vector  space  defines  a convex, 
compact  set  with  symmetry  in  each  of  the  miss  vector  component 
directions.  The  miss  vector  set  so  established  can  be  related  to  the 
statistical  expectation  of  the  various  system  error  sources.  The  INS 
allocation  of  this  acceptable  expected  error  set  can  be  reflected  onto 
the  INS  error  state  vector  coordinates.  Since  all  INS  error  state 
coordinates  may  not  contribute  directly  to  system  miss,  the  INS 
error  state  set  so  defined  may  not  be  bounded.  With  our  previous 
linear  assumption  it  will  be  convex  and  symmetric  with  respect  to  all 
component  directions.  Thus,  we  have  a set  which  we  can  view  as  a 
restriction  on  the  INS  error  state  covariance  matrix. 

We  shall  not  treat  system  cost  functions  or  miss  vector  coordi- 
nates in  further  analyses  since,  under  our  assumptions,  these  concerns 
can  be  stated  as  constraints  on  the  INS  error  state  covariance  matrix. 
The  INS  error  state  covariance  restrictions  so  derived  are,  in  turn, 
distributed  across  the  INS  component  error  sources --the  gravitational 
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model  being  one  such  element.  In  this  manner,  we  view  the  design 
criteria  for  the  gravitational  model  to  be  stated  in  terms  of  bounds  on 
the  INS  error  state  covariance  matrix  resulting  from  gravitational 
modeling  errors  acting  alone.  An  inherent  assumption  in  this  process 
is  the  independence  of  gravitational  modeling  errors  from  other  INS 
error  sources  (e.g.  accelerometer  noise)  and  from  other  system 
error  sources  (e.g.  errors  in  the  estimate  of  ballistic  coefficient). 

With  this  manner  of  specifying  the  system  level  performance 
of  the  gravitational  model,  we  need  to  develop  analytical  means  of 
estimating  the  INS  error  covariance  matrix  which  results  from  gravi- 
tational modeling  errors  acting  alone.  The  necessary  elements  of 
this  analysis  have  been  discussed  previously.  We  need  the  mission 
type  description  to  fix  the  INS  error  propagation  model.  The  geo- 
graphic region  defined  by  the  mission  type  can  be  used  along  with 
gravity  survey  data  to  estimate  the  covariance  of  anomalous  gravita- 
tion. We  must  devise  a means  of  altering  this  anomaly  covariance  or 
anomalous  potential  covariance  to  include  the  effects  of  our  additional 
modeling  of  the  gravitational  disturbance.  That  is,  we  need  to  derive 
the  covariance  of  the  residual  field  qu  itity. 

Then,  an  analytical  means  must  be  derived  to  propagate  this 
expected  gravitational  modeling  residual  noise  through  the  INS  taking 
into  account  the  design  mission  effects  on  both  the  residual  gravitation 
field  statistics  and  on  the  INS  error  propagation.  The  means  for 
accomplishing  this  task  must  be  statistical.  If  we  knew  the  gravitational 
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field  exactly  along  all  mission  trajectories,  we  could  form  the 
ensemble  statistics  through  well-known  statistical  techniques.  Lack- 
ing such  extensive  information,  we  can  form  the  outer  product  of  the 
INS  error  vector  and  apply  our  expectation  operation  over  the  mission 
region.  As  discussed  in  Appendix  C,  the  INS  error  covariance  from 
gravitational  modeling  errors  has  been  studied  by  simulation  methods 
for  some  rather  benign  INS  missions.  Since  we  want  a method  which 
is  generally  applicable,  we  need  analytical  methods  which  are  not 
restricted  to  simple  mission  geometries. 

This  generalized  design  mission  analysis  technique  will  pro- 
vide the  answer  to  whether  a proposed  model  concept  and  complexity 
level  provide  acceptable  system  level  performance.  The  system 
resource  costs  can  then  be  compared  for  various  model  concepts  or 
the  system  level  performance  can  be  compared  for  fixed  system 
resource  investmmt.  The  design  mission  context  allows  us  to  address 
the  overall  system  mission  performance  by  analyzing  a finite,  but 
representative,  subset  of  missions.  The  information  on  system  per- 
formance and  system  resource  costs  can  thus  be  made  available  to 
answer  the  design  trade-off  questions  which  will  arise  in  selecting  an 
extended  gravitational  model. 

B.  Specific  Objective 

The  goal  of  this  research  is  to  provide  analytical  methods  for 
selecting  a computationally  efficient,  INS  gravitational  model  which  is 
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accurate  in  terms  of  expected  INS  accuracy.  The  research  will  focus 
on  providing  a means  of  estimating  INS  errors  due  to  gravitational 
model  errors  acting  alone.  The  system  resource  cost  will  be  approx- 
imated as  a linear  function  of  the  number  of  gravitational  model 
parameters. 

C.  Implementation  Study 

The  method  discussed  above  will  vary  with  the  model  concept. 
The  reason  is  that  the  residual  field  depends  entirely  on  the  manner 
in  which  the  anomalous  field  is  modeled.  The  Poisson  Integral  (36) 
will  be  analyzed  as  an  example  to  demonstrate  the  application  of  the 
method.  The  hypothetical  design  study  will  be  to  determine  the 
Poisson  Integral  approximation  grid  which  meets  an  INS  accuracy 
criteria  with  the  minimum  number  of  grid  elements  (i.  e.  number  of 
parameters).  The  Poisson  Integral  was  selected  as  an  example 
because  (1)  it  has  not  received  much  attention  in  previous  works;  (E) 
current  satellite  altimetry  will  provide  the  necessary  data  base  over 
large  regions  of  the  Earth;  (3)  this  integral  may  prove  amenable  to 
solution  by  transform  techniques;  and  (4)  the  relationship  between 
gravitation  and  potential  is  more  straightforward  than  the  relationship 
to  gravitational  anomaly  which  is  the  data  basis  for  Stokes  Integral. 

D.  Assumptions 

The  following  assumptions  will  be  made  in  this  research.  The 
assumptions  fall  into  one  of  two  categories: 


74 


1.  Those  assumptions  which  provide  details  outside  the  direct 
path  of  the  proposed  research  (e.g.  on  survey  errors),  and, 

2.  Those  assumptions  which  clarify  the  objective  and  ensure 
mathematically  tractable  results. 

The  proposed  research  will  not  deal  directly  with  either  survey 
errors  or  total  gravitational  model  algorithm  computational  require- 
ments. These  quantities  are  needed  to  complete  the  analysis  and  we 
shall  assume 

1.  That  system  resource  cost  will  be  reflected  in  the  number  of 
model  parameters  (i. e.  number  of  grid  elements  for  Poisson 
Integral  example), 

2.  That  gravity  survey  and  subsequent  data  processing  will 
supply  anomalous  potential  estimates  on  an  arbitrary,  compu- 
tationally convenient  grid,  and 

3.  That  survey  errors  are  independent  of  the  data  and  are  a 
white,  Gaussian,  two-dimensional  sequence. 

The  assumptions  which  will  be  made  to  yield  tractable  mathe- 
matics have,  for  the  most  part,  been  discussed  previously.  The 
additional  assumptions  on  the  anomalous  and  residual  fields  seem 
reasonable  and  yield  significant  simplifications  in  the  covariance 
functions  {see  Appendix  C).  These  assumptions  are 

1.  That  system  performance  cost  will  be  stated  as  a restric- 
tion on  the  INS  error  state  covariance  matrix  due  to  gravitational 
model  errors  acting  alone. 
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2.  That  the  design  mission  represents  the  entire  mission 
region  for  INS  error  propagation  modeling, 

3.  That  the  statistics  of  the  anomalous  gravitational  field  are 
known  or  can  be  approximated  over  the  reference  surface  in  the  geo- 
graphic area  of  the  mission  region,  and 

4.  That  the  statistics  of  the  anomalous  or  residual  field  are 
well-approximated  by  assuming  either  field  is  a homogeneous  and 
isotropic  process  over  the  reference  surface. 

E.  Proposed  Approach 

A major  portion  of  the  approach  to  this  problem  is  contained  in 
the  Problem  Context  section.  For  a general  class  of  INS  applications, 
the  overall  model  selection  process  is  based  on  the  performance 
results  from  a finite  set  of  design  missions.  The  basic  problem 
remaining  is  to  develop  an  analytical  procedure  for  determining  the 
expected  system  performance,  on  any  specific  design  mission.  To  us 
system  performance  is  synonymous  with  the  INS  error  state  covari- 
ance matrix,  given  (1)  the  model  concept,  (2)  the  model  complexity 
level,  (3)  the  regional  anomalous  field  statistics,  (4)  the  design  mis- 
sion, and  (5)  the  survey  error  statistics. 

This  problem  will  be  broken  into  two  distinct  parts  (see  Figure 

5).  In  the  first  phase,  we  shall'characterize  the  residual  gravitational 

field  statistically.  For  our  Poisson  Integral  problem  this  summary 

will  be  in  the  form  of  the  residual  potential  covariance  function  over 

the  reference  surface.  This  residual  field  statistic  will  contain  all 
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necessary  information  on  the  gravitational  model,  the  mission  region 
and  the  survey  errors.  The  second  phase  of  this  analysis  is  generally 
applicable  for  all  model  concepts.  Given  the  residual  field  statistic 
from  the  first  phase,  we  shall  form  a statistical  estimate  of  the  system 
cost  function  based  on  the  design  mission  and  our  INS  error  propaga- 
tion model  (Refs  1 or  2). 

Along  with  this  system  performance  cost  our  implementation 
study  will  use  the  model  complexity  as  the  system  resource  cost.  We 
shall  investigate  the  solution  to  the  design  trade-off  question:  Find 
the  least  complex  Poisson  Integral  approximation  which  meets  an  INS 
accuracy  level  specification. 

The  start  of  this  overall  analysis  is  to  form  the  residual  field 
statistics.  For  our  Poisson  Integral  problem,  we  assume  that  the 
anomalous  potential  covariance  function  K(^,jR’)  is  known  (see 
Appendix  C for  a general  discussion  of  such  covariance  functions). 

The  radius  vectors  ^ and  are  on  the  reference  surface  of  radius  R. 
The  covariance  argument  will  usually  be  simply  the  central  angle, 
separating  ^ and  the  form  K(^,^*)  will  be  used,  however,  to  allow 
a more  general  class  of  arguments.  We  shall  use  K(^,^')  for  the 
worldwide  covariance  function.  The  symbol  0 will  signify  the  mission 
region  being  considered  and  K(^,^';0)  will  be  the  anomalous  potential 
covariance  formed  from  data  in  the  geographic  region  associated 
with  0 . 
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The  modeling  of  anomalous  gravitation  results  in  a residual 
field  with  statistics  dependent  on  the  anomalous  field  statistics  and  the 
model  approximation  performance.  The  residual  potential  covariance 

A 

function  will  be  called  K (^,  ^')  and  this  quantity  will  be  formed  for 
our  Poisson  Integral  example. 

We  are  using  this  example  to  provide  concrete  terms  for  what 
would  otherwise  be  an  abstract  analysis.  The  nature  of  how  we 
approximate  the  Poisson  Integral  will  influence  the  results.  Recall 
(36) 
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where  T is  shown  in  Figure  4 on  page  56.  We  shall  approximate  (36) 
directly  rather  than  use  the  Equation  (44).  This  choice  was  made  for 
example  simplicity  since  otherwise  we  would  need  to  analyze  the  three 
integrals  of  (44)  in  addition  to  (36).  The  spatial  derivative  will  be 
approximated  by  a discrete  number  of  evaluations  of  (51)  in  the  vicinity 
of  the  navigation  position  estimate,  £.  We  shall  approximate  (51)  with 
a finite  sum  replacing  the  integral.  The  sum  will  be  formed  from  a 
uniform  grid  of  square  surface  elements.  We  shall  neglect  any  prob- 
lems these  square  elements  encounter  in  covering  the  spherical  refer- 
ence surface  because  the  grid  extent  will  be  limited.  From  our 
previous  Stokes  Integral  discussion,  we  expect  the  integral  approxi- 
mation to  converge  slowly  as  the  grid  extent  increase  past  ten  degrees 
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of  central  angle  from  the  evaluation  point.  Depending  on  the  stringency 
of  our  INS  accuracy  requirement,  we  can  anticipate  a relatively  small 
grid  will  suffice  for  many  applications.  The  approximation  to  (51) 
will  be 


T(lO  = ) Y* 

4n 

3=1 


IRj-r|3 


(52) 


where  A*  is  the  area  of  the  surface  element  and  m is  the  number 
of  grid  elements.  With  our  uniform  grid  we  get 


^(rj  . R(r^-R^)  A^  t{R:) 
4jt  ^ — — 


3=1  —3  - 


(53) 


The  Poisson  Integral  approximation  model  concept  has  com- 
plexity level  indicated  by  m in  (53) --the  number  of  grid  elements. 
Since  the  radius  of  the  reference  surface  is  fixed  at  R,  we  need  two 
parameters  to  specify  the  grid  point  which  applies  to  the  survey 
data  ^ (Rj).  Thus,  three  computer  memory  storage  locations  might 
be  required  to  identify  the  required  information.  We  can  use  the  uni- 
form grid  structure  to  reduce  this  storage  to  just  the  anomalous 
potential  data  with  a grid  center  point,  for  example. 

The  kernel  of  (51)  is  well  behaved  for  all  ^ such  that  r > R.  So 

if  our  data  | T (^i)|.  is  without  survey  error,  we  expect  (52)  or  (53) 
J j = l 

to  converge  to  the  true  value  in  the  limit  as  ni  and  each  Aj  -^O. 
This  quality  is  the  prime  motivation  for  choosing  such  integrals  as 
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the  basis  for  a model  concept.  Another  way  of  specifying  the  grid  is 
by  the  spherical  cap  size  and  the  grii  element  dimension. 

The  spherical  cap  is  the  reference  surface  area  within  a central 
angle  (discussed  previously)  of  the  evaluation  point  £.  The  element, 
being  square,  is  specified  by  the  length  of  the  side,  h.  We  expect 
convergence  as  ^ and  as  h ->0.  The  choice  of  this  pair  is  moti- 

vated by  the  symmetry  of  the  (51)  kernel  radially  from  jr  and  by  the 
connection  of  h with  the  Shannon  Sampling  Theorem- -here  applied  to 
representation  rather  than  sampling. 

Clearly  for  a fixed  m,  there  exists  a family  of  ('/'q,  h)  pairs 
which  satisfy  the  equation 


m 


total  cap  area 
grid  element  area 


2 R^(l-cos  >j/o) 


(54) 


We  envision  the  INS  specification  as  a partition  on  the  plane 

dividing  the  possible  grid  designs  into  acceptable  and  not  acceptable 
grid  specifications.  Our  task  in  the  implementation  study  will  be  to 
select  from  the  "acceptable"  set  the  grid  specifications  which  mini- 
mize m. 

We  shall  discuss  this  problem  later,  our  task  now  is  to  form  the 
K function  and  this  grid  design  data  is  needed.  Recall  that  in  forming 
K*  we  assume  the  grid  design  (i,  e.  and  h)  is  given  as  part  of  the 
model  definition.  Clearly  the  statistic  K is  a key  to  our  predicting 
INS  error  performance  with  our  compensated  gravitational  model. 
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The  method  of  determining  from  K(I^^';0)  nnd  the  grid 

specifications  is  a major  problem  to  be  undertaken.  We  know,  by 
Shannon's  Sampling  Theorem,  that  the  grid  spacing  affects  the  fre- 
quency content  of  the  residual  field.  So,  expanding  K(^, ^';0)  in 
spherical  harmonics  may  provide  the  insight  for  accounting  for  grid 
size  effects.  The  cap  size  study  should  follow  the  previously  dis- 
cussed methods  of  Hirvonen  and  Moritz  (Ref  64  and  Ref  14:242-243). 

Since  the  survey  errors  are  independent  from  the  data  (potential 

in  our  Poisson  Integral  example),  we  can  treat  their  contribution  at 

- ^ 
the  covariance  matrix  level  and  write  K as 

K*(R.R';0)  = K^(R,  R';0)  + Kg(R,  R-)  (55) 

where  is  the  covariance  of  the  residual  field  assuming  perfect 
data  and  Kg(_R,JR')  is  the  extended  survey  error  covariance  function 
which  results  from  the  Poisson  Integral  acting  on  Ks(^i'^j)* 

We  want  Kg  and  Kj^  to  be  symmetric  with  respect  to  the  central 
angle  argument  for  later  computational  ease.  We  can  study  Kg  by 
considering  the  Poisson  Integral  as  merely  performing  an  interpolation 
between  grid  points.  The  details  for  this  development  have  not  been 
completed.  A simple  one  dimensional  example  will  illustrate  the 
intended  approach.  Suppose  our  grid  points  are  Rj  and  R£  and  we  want 
an  evaluation  at  Re[R^,R2].  Then  we  may  express  R as 

R = aRj  + (1  - O')  R2  (56) 
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for  0 5 ® 5 1,  Using  the  parameter  a and  the  grid  data  we  form  our 
interpolation  as 


•^(R)  =a^(Ri)  + (l  -o)^(R2) 


(57) 


Let 


T(Ri)  = T{Ri)  + qi 


(58) 


for  i = 1 and  2 where  is  the  uncorrelated  error  associated  with  the 
survey  data  point.  We  have  then 


( 2 , . . 
(Tg  for  1 = j 

0 else 


(59) 


Using  our  definition  for  error  again,  let 

^‘(R)  = T(R)  + q^(R)  (60) 

where  q.j.  represents  the  total  error  at  R.  Combining  (57),  (58)  and 
(60)  we  get  q^(R)  = 9jj^(R)  + qs(^)  (^>1) 

where  qjjj(J^)  is  the  error  due  to  our  interpolative  model  given  by 


q^(R)  =«T{Rj)  + (1  -a)T(R?j  - T(R)  (62) 

and  where  qg(R)  is  the  error  at  R due  to  survey  errors  at  Rj  and  R2 
given  by 

qg(R)  = «qi  + (1  -«)  q2*  (53) 
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. u 0>m9^ 


By  our  assumed  independence  of  the  survey  errors  with  the  true  field 
we  get 

t + ^[q^R)]*  (64) 

Now  the  survey  error  contribution  to  this  covariance  can  be  calcu- 
lated using  the  property  expressed  in  (59): 

^[q^R)]  .=  (l  - 2c^+2ff^)(rg^.  (65) 

With  <r“  a constant  statistic  across  the  grid,  we  would  find  this  form 
repeated  between  each  successive  pair  of  grid  points.  The  Poisson 
Integral  approximation  can  be  viewed  as  a two  dimensional  interpola- 
tion using  all  grid  data  points.  While  the  analysis  for  Kg(^,jR'f  of  (55) 
will  be  more  complicated,  the  example  above  demonstrates  the  intended 
approach. 

The  statistics  for  9jn(R)  above  would  be  ^[q^(R)].  This  quantity 

corresponds  to  K^(^^';0)  of  (“5).  We  shall  approach  this  quantity 

in  an  entirely  different  man  .''r.  The  grid  direction  will  be  such  that 

the  downrange  groundtrack  lies  along  a partition  boundary.  The 

distance  between  data  points  will  be  h so  we  can  apply  sampling  and 

approximation  theory  to  form  K^(^,^';©)  from  K{R,^';0).  This 

analysis  is  incomplete  at  this  time,  but  the  approach  will  be  to  assume 

Km(^i»^i'®)  “ ® since  the  infinite  kernel  of  (5Z)  or  (53)  implies  the 
J J 

weighting  cf  the  true  field  at  that  point  will  be  such  that  the  remainder 
of  the  data  points  can  be  ignored  (recall  deals  with  perfect  data 
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since  incorporates  the  survey  errors).  A means  of  formulating 
the  error  at  points  on  the  reference  surface  between  the  gridded  data 
points  is  required.  Given  this  error  function,  the  covariance  will  be 
given  by  taking  an  expectation  over  all  gravitational  fields  we  expect 
to  encounter.  This  expectation  can  be  construed  to  be  over  © and 
K(R,^';©)  should  provide  all  the  required  information. 

Assumcng  Kg  and  can  be  calculated  or  approximated,  we 
have  K by  (55)  which  summarizes  the  statistics  of  the  residual  poten- 
tial field  throughout  0.  Let  9 signify  a mission  trajectory;  then  our 
analysis  must  consider  all  0€©.  Let  6^  signify  the  design  mission 
which  represents  the  mission  region  0.  We  can  now  form  an  estimate 
of  the  INS  error  state  covariance  matrix. 

Using  Reference  1 methods,  we  can  model  the  INS  error  propaga- 
tion for  our  gravitational  residual  as 

X(t;fl)  = F(t;<?)X(t;e)  + G(t;ff)  u(t;6)  (66) 

where  X is  the  INS  error  state  vector  including  position,  velocity,  and 
attitude  errors;  F is  the  error  propagation  matrix  dependent  on  the 
time  in  to  trajectory  6 by  the  associated  geography  and  geometry;  G is 
the  driving  noise  distribution  matrix;  and  u is  the  particular  residual 
field  errors  along  9 . The  driving  noise  £will  contain  the  residual  ' 
acceleration  ^ and  possibly  the  residual  potential  ST  for  a barometric 
altimeter  aided  INS,  For  completeness  we  assume 


84 


u = 


S 

8T 


(67) 


Since  our  analysis  is  for  gravitational  modeling  errors  only,  we 
assume  initial  conditions  for  (66)  of  }C(O;0)  = 0.  (68) 

We  can  form  a solution  to  (66)  using  a state  transition  matrix 
which  satisfies 

i(t,ti;fl)  = F(t;fl)<I)(t,tj;e)  , (69) 

and 

<l)(t,t;fl)=I  (70) 

where  I is  the  identity  matrix.  Such  a matrix  also  satisfies  a semi- 
group property 


:=<t>(t3,t2;6)<l>(t2,ti;e) 

Using  <t(t,  tjjff)  with  (66)  and  (68)  we  get, 

.t 


(71) 


X(t;9)  =y*4)(t,p;5)G(p;e)u(p;0)dp 


(72) 


Forming  the  outer  product  of  this  error  vector  yields 

.t.t 


X(t;fl)x"^(t;9)  = jJ"  <J)  (t,p;fl)G(p;0  )u(p;e)u'^(q;9)G'^(q;0  )<l>'^(t,q;0)dpdq 

(73) 


o o 


Now,  we  define  the  INS  error  state  covariance  function,  as  the 

expectation  over  all  missions  9 within  the  mission  region  0. 
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That  is, 

Pvx^t)  = e [X(t;fl)X'^(t;fl)].  (74) 

9eQ 

Equations  (74)  and  (73)  give  a truly  general  definition  of  Pj^.  due  to 
gravitational  model  errors.  For  cases  with  a small  set  of  trajectories 
and  with  a well  known  gravitational  field,  these  equations  provide  the 
total  analysis.  For  most  cases  of  interest  the  mission  region  0 con- 
tains a continuum  of  trajectories  and  we  know  very  little  about  the 
field  for  any  particular  fl  . We  seek  then  a less  general  form  of  (74) 
which  is  better  aligned  with  the  information  we  expect  to  have  avail- 
able. 

One  mission  . pace  partition  requirement  was  that  each  mission 
region,  such  as  0 being  discussed  here,  is  well  represented  by  a 
single  design  mission  9^  as  far  as  INS  error  dynamics  are  concerned. 
With  this  stipulation  applied  to  (73)  we  get 

X{ti9)^{t;9)  ~JJ <t(t,p)G(p)u(p;^;)u'^(q;  0)G’^(q)<l>'^(t,q)dpdq  (75) 

o o 


where  ^ and  G are  defined  for  6^.  Now  applying  (74)  to  (75)  we  can 
take  the  expectation  operation  inside  the  integral  to  yield 


where 
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(77) 


Q[r(p),r(q),K*(R»R’;e)]  = e |u(p;0)u'^(p;0)| 

dee 

and  £(p)  and  £(q)  are  on  Oq  since  0^  also  represents  the  geometry 
associated  with  the  OeQ.  The  Q- matrix  is  discussed  in  Appendix  E. 
The  covariance  function  K {R,R';Q)  characterizes  the  residual  field 
statistics  as  a function  of  the  geometry  inherent  in  the  ^ and  R*  argu- 
ments. This  function  can  be  continued  upward  (see  Appendix  C)  to 
form  K(_r,£';6).  Now  given  3r(p;  6)  and  r^(q,  5)  as  geometric  arguments 
associated  with  the  trajectory  $,  we  can  form  K [£(p, ^ ),£(q, ^ );  ©]. 
From  this  fvinction  we  can  form  a Q-matrix  for  each  6.  Since  the  K* 
func  .ion  is  a function  of  the  geometry  of  and'£(q,  0)  we  can  use 

the  evaluation  along  the  design  mission  0^  to  approximate  the  expecta- 
tion over  all  0€9. 

In  (76)  we  have  the  analysis  method  we  have  been  seeking.  The 
covariance  provides  the  data  for  either  computing  system-level 
performance  cost  (see  Appendix  D)  or  for  comparison  to  an  INS-level 
accuracy  criteria  for  the  gravitational  model.  The  gravitational 
modeling  errors  over0  are  summarized  in  the  Q-matrix  and  their 
expected  effect  on  INS  estimates  are  propagated  through  the  G and  <t> 
functions.  Equation  (76)  appears  to  be  a useful  form  for  other 
analyses,  as  well.  For  example,  suppose  our  original  question  were 
"Do  we  need  to  consider  more  than  the  reference  field  model  to  meet 
our  INS  gravitational  model  error  budget?"  To  answer  this  question 
we  apply  (76)  to  each  design  mission  using  K(^,^';0)  rather  than  K* 
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since  our  hypothesis  is  only  the  reference  field.  We  can  estimate  the 
ultimate  accuracy  attainable  by  forming  an  equivalent  continuum 
version  of  Kg  and  applying  (76)  with  this  covariance  under  an  assumed 
perfec*:  Poisson  Integral  approximation.  These  analyses  form  a worst 
case,  with  0 ),  and  a best  case,  with  Kg(^,^'),  set  of  numbers 

which  will  be  useful  in  judging  the  relative  effectiveness  of  our  model. 

Also,  if  one  were  designing  a Kalman  filter  for  an  INS  applica- 
tion, (76)  allows  for  explicit  consideration  of  this  system  noise  source. 
Another  possible  use,  along  this  line,  is  to  let  (76)  provide  the  sta- 
tistics for  a truth  model  gravitational  noise  model  at  the  INS  estimate 
level. 


Equations  of  the  form  of  (76)  appear  in  Kalman  filter  develop- 
ment and  are  typically  put  in  differential  equation  form  for  computer 
solution.  The  appearance  of  r(p)  and  £(q)  in  the  Q-matrix  argument 
list  will  complicate  this  approach.  The  first  step  is  to  differentiate 
(76)  with  respect  to  the  mission  time  t: 


Pj^(t)  =G(t) y*  Q[r(t),r(q);K*]G'^(q)<X>'^(t,q)dq 

+ |y’^<l>(t,p)G(p)Q[r(p),r(t);K*]dp|GT(t) 


+ F(t)Pj^^(t)  + P^(t)pT(t) 


(78) 


subject  to 

Pxx^O)  = [0]-  * (79) 
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Q is  a covariance  matrix  hence  symmetric  so 


Q[r(t),r(q);K*l  =Q[r(q),r(t);K*l 


(80) 


We  can  rewrite  (78)  as 


= G(t)D(t)  + DT(t)cFt)  + Pxx(t)F'^(t)  + F(t)Pxx(t)  (81) 


where 


D(t)  = y*Q[r(t),r{q);K*]GT(q)<l>'^(t,q)dq 


(82) 


The  £(t)  in  the  Q argument  of  (82)  makes  it  undesirable  to  continue 
with  a D(t)  equation  for  computational  purposes.  The  set  (81),  (82) 
and  (69)  subject  to  initial  conditions  of  (70)  and  (79)  form  a complete 
set  for  computational  analysis.  The  approach  will  be  to  write  com- 
puter programs  to  solve  the  differential- integral  equation  set  so 
defined  as  a means  of  estimating  the  INS-level  performance  of  the 
gravitational  model. 

While  the  differential  form  of  (82)  is  computationally  unwieldy, 
this  equation  is  importaut  for  the  insights  it  affords. 


i)(t)  = Q[r(t),  r(t);K*]GT(t)  + D(t)F^(t) 

+ y j~Q[r.,£(q):K*]!  v(t)G'^(q)(J>'^(t,q)dq. 


(83) 


The  first  term  on  the  right  hand  side  of  (83)  represents  the  new 
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gravitational  disturbance  entering  at  point  £(t).  The  second  term  is 
just  the  propagation  of  previous  errors.  The  third  term  represents 
a "missing  link"  between  the  analysis  leading  to  (76)  and  our  previous 
intuitive  examples.  The  partial  of  Q~with  respect  to  a spatial  variable 
can  be  shown  to  be  related  to  the  gravitational  tensor  F of  the  residual 
field.  The  v(t)  term  is  the  mission  velocity  which  in  conjunction  with 
informs  the  time  varying  gravitational  disturbance  derivative. 

In  summary,  the  system-level  performance  cost  estimation 
occurs  in  three  distinct  steps.  This  process  is  shown  in  flow  chart 
form  in  Figure  5.  The  model  concept,  mission  region,  and  survey 
characteristics  are  used  to  form  a residual  field  covariance  function. 
This  function  along  with  the  INS  error  dynamics  and  the  design  mission 
are  used  to  form  the  INS  error  state  covariance  matrix  as  a fmiction 
of  time  into  the  mission.  This  covariance  function,  contains 

all  the  information  needed  to  calculate  the  expected  system  cost  (see 
Appendix  D);  recall  that  we  do  not  plan  to  use  J(t)  in  our  example  but 
shall  assume  the  requirements  placed  on  ^[J(t)]  be  interpreted  as 
restrictions  on 

With  this  presentation  of  the  approach  for  performance  cost 
evaluation,  let  us  turn  to  the  implementation  study.  V/e  have  (76)  and 
(541  to  represent  the  opposing  performance  and  resource  costs.  Our 
example  application  is  to  find  the  minimal  grid  (least  m)  which  meets 
an  INS  accuracy  requirement.  This  performance  requirement  might 
be  to  arrive  at  terminal  time  t|  with  an  rms  position  uncertainty  due 


Figure  5.  Gravitational  Modeling  Accuracy  Cost 
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to  the  gravitational  model  of  no  more  than  Z feet. 

The  control  parameters  for  grid  specification  are  and  h. 

Our  task  is  to  select  and  h to  minimize  m from  (54)  such  that  the 
performance  cost  from  (76)  meets  specification.  This  problem  falls 
into  the  category  of  constrained  optimization.  We  know  performance 
■will  improve  as  is  increased  since  more  information  is  available. 
We  expect  performance  to  improve  as  h decreases,  however  the  inter- 
play between  h,  the  mission  velocity  profile,  and  the  Schuler  loop  may 
cause  a surprise.  With  either  decreasing  h or  increasing  we  get  an 
increase  in  m-system  resource  cost.  Therefore,  we  can  expect  the 
minimal  grid  solution  to  lie  on  the  boundary  of  the  hypothetical  per- 
formance cost  partition  on  the  '#'Q-h  plane. 

The  example  search  logic  is  shown  in  Figure  6.  This  figure 
is  designed  to  form  an  overall  flow  chart  with  Figure  5.  The  actual 
search  logic  will  be  defined  after  some  characterization  studies;  e.g. 
constant  performance  cost  and  constant  resource  cost  studies.  These 
studies  should  provide  information  on  the  smoothness  of  constant  per- 
formance curves  on  the  '/'g-h  plane.  We  can  also  gain  some  insight 
by  studying  the  variations  in  performance  cost  along  a constant  m 
line. 

This  concludes  the  description  of  the  Proposed  Approach.  Any 
analysis  must  face  a test  to  prove  its  worth.  With  this  in  mind,  we 
direct  the  discussion  to  the  methods  which  will  demonstrate  the 
effectiveness  of  these  cinalyses. 
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Figure  6,  Minimal  Grid  Search  Logic 


F.  Confirmation  Method 


The  manner  of  evaluating  the  developed  theory  -will  be  by  simu- 
lating true  gravitational  fields,  by  applying  the  theory  to  these  fields 
for  a predefined  mission,  and  by  studying  the  errors  in  a simulated 
INS  using  the  derived  compensation  model.  Two  types  of  simulated 
gravitational  fields  come  to  mind.  One  type  is  pure  simulation.  An 
analytical  field  is  constructed  to  reflect  or  approximate  the  statistics 
of  the  observed  field--say  by  producing  the  same  K(^,  R';0)  as  an 
empirically  derived  function.  Another  simulation  approach  is  to 
attempt  to  reproduce  the  field  in  a general  area  by  closed  form  approx- 
imations (e.g.  point  masses  as  in  Ref  24).  Both  methods  have  advant- 
ages. The  pure  simulation  allows  us  to  control  the  field  to  the  (nearly) 
identical  statistics  upon  which  our  analysis  will  be  based.  The  real 
field  approximation  will  allow  a more  convincing  test  of  our  hypothesis. 

To  use  the  analytical  fields  the  survey  specifications  will  be 
solved  for  a predetermined  mission  type.  Grid  sampling  will  be  accom 
plished  by  direct  computation  of  the  analytical  field  equations.  Survey 
errors  can  be  added  to  these  data  for  Monte  Carlo  simulations.  Since 
the  gravitational  modeling  errors  are  the  issue,  the  simulated  INS  will 
be  initialized  with  perfect  knowledge  of  the  initial  state  vector  and  with 
perfect  alignment.  Simulated  mission  position  and  velocity  along  with 
the  simulated  true  gravitational  field  will  be  used  to  derive  INS  sensor 
inputs.  The  INS  algorithm,  not  just  the  error  propagation  model,  will 
be  programmed  with  the  INS  reference  gravitation  field  compensated 
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with  the  Poisson  Integral  approximation  of  gravitational  disturbance. 

The  resulting  INS  errors  can  be  subducted  from  the  simulated 
true  and  the  INS  estimated  states.  These  error  time  histories  can  be 
compared  with  our  specification.  A Monte  Carlo  study  with  variations 
on  the  true  field  and  the  trajectory  fixed  or  variations  on  the  trajectory 
with  the  field  fixed  is  necessary  to  demonstrate  the  required  ensemble 
performance. 

G.  Ancillary  Questions 

The  confirmation  methods  will  measure  our  ability  to  construct 
grids  which  deliver  the  advertised  INS  accuracy.  The  Poisson  Integral 
will  require  a major  investment  of  our  computer  resources,  and  it  is 
reasonable  to  ask  what  we  are  getting  in  return.  For  example,  can 
less  sophisticated,  but  perhaps  less  flexible,  compensations  be  com- 
pared. To  put  the  Poisson  Integral  in  perspective  this  question  and 
others  will  be  addressed: 

1.  How  effective  is  a pre-mission  initial  condition  offset  based 
on  the  nominal  trajectory  and  a possible  more  complex 
model  for  ground-based  evaluation? 

2.  What  is  the  sensitivity  of  our  performance  to  the  choice  of 
nominal  trajectory? 

3.  How  sensitive  are  our  results  to  errors  in  covariance  func- 
tion K(R,  R')? 

4.  Can  the  integral  be  implemented  by  transform  techniques? 
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• IV . Outline  of  the  Final  Report 

The  dissertation  will  be  based  on  this  prospectus.  The  back- 
ground provided  herein  will,  with  some  modification,  be  repeated  in 
the  introductory  discussion.  The  completion  of  the  theoretical  analysis 
with  the  necessary  basis  from  the  proposed  approach  herein  will  be 
recorded.  The  numerical  confirmation  tests  and  ancillary  questions 
will  be  provided  in  the  results  section.  The  major  subdivisions  and 
their  general  content  should  be  as  follows. 

1.  Introduction;  The  background  discussion  necessary  to  put 
the  theoretical  work  in  context.  The  reference  gravity  model 
will  be  more  thoroughly  discussed--perhaps  as  an  appendix. 

2.  Theory  of  System- Level  Performance  Cost  Estimation  for 
Gravitational  Models;  Completion  of  the  development  leading 
to  (76). 

Theory  of  Poisson  Integral  Approximation  (for  INS  applica- 
tion);  The  fundamental  theory  will  be  cited  and  summarized 
for  the  Poisson  Integral  and  for  the  approximation  of  such 
integrals  by  finite  element  summation.  The  developments  of 
the  Proposed  Approach  will  be  repeated  and  completed. 

4,  Implementation;  Those  considerations  necessary  to  bridge 
from  theoretical  development  to  computer  code. 

5.  Example  Application  Problem;  Discussion  of  the  minimal 
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grid  problem. 

6.  Numerical  Results;  The  confirmation  studies  and  replies 
to  the  ancillary  questions. 

7.  Summary  and  Conclusions;  Self-explanatory. 

8.  Suggestions  for  Future  Research;  Self-explanatory. 

9.  Appendices;  As  required  to  package  lengthy  developments 
not  directly  in  line  with  the  theoretical  development. 
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V.  Air  Force  Applicability 

The  applicability  of  the  proposed  research  to  the  Air  Force  has 
been  touched  on  in  several  sections  (e.g.  pages  24-25).  The  precise 
amd  unaided  INS  is  an  obvious  flexible  asset  in  employing  strategic 
weapons.  The  techniques  proposed  herein,  once  developed,  could  also 
be  applied  to  the  test  flight  problem.  The  regression  analyses  for 
identifying  INS  component  error  model  parameters  would  employ  the 
gravitational  compensation  model  to  remove  the  corruption  of  the 
results  due  to  the  spectral  spillover  of  gravitati  mal  disturbances. 

The  potential  for  this  theory  being  directly  applied  cannot  truly 
be  estimated.  If  the  upward  continuation  integrals  can  be  cast  in  a 
form  compatible  with  Fast  Fourier  Transform  (FFT),  the  potential 
application  increases.  The  onboard  memory  space  in  most  operational 
systems  will  make  the  expensive  fine  model  unacceptable.  With  FFT 
formulation  we  can  consider  an  off-line,  separate  gravitation  com- 
puter (see  Figure  7)  which  is  an  input-output  device  on  the  data  bus. 
Hardware  is  presently  available  to  perform  the  FFT  "butterfly"  opera- 
tions. This  concept  would  be  an  asset  for  strategic  bomber  and 
missile  application.  Even  in  tactical  situations  we  might  find  a highly 
accurate,  all-inertial  capability  useful  in  many  ways. 
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Figure  7.  Future  Application  Possibility 


VI.  Originality 


The  proposed  approach  is  original  in  the  concept  of  focusing  on 
the  INS  results  rather  than  anomalous  gravitation  per  se.  The  idea 
of  using  an  INS  error  budget  as  a criteria  or  design  constraint  is 
common  practice,  but  the' application  to  this  problem  of  gravitational 
modeling  is  new.  Equation  (76)  and  the  rationale  leading  up  to  it  are 
also  original.  The  general  applicability  of  (81),  (82)  and  (69)  mean 
they  can  be  employed  in  many  ways  whenever  the  rationale  fits  and 
the  covariance  function  is  known  or  approximated.  Applying  these  con- 
cepts to  the  grid  specification  problem  is  an  uncompleted  task.  The 
end  result  will  be  a definite  extension  of  grid  specification  which  might 
have  application  in  survey  design  tasks  as  well.  The  remaining  task  of 
deriving  the  residual  potential  covariance,  K {^,R*;Q)  from  the  anoma- 
lous potential  covariance,  K(^,  R';0),and  the  grid  specifications  should 
also  provide  an  original  contribution. 


Appendix  A 
The  Gravity  Tensor 

The  gravity  tensor  plays  an  important  role  in  the  INS  error 
propagation  model,  and  it  has  received  much  attention  recently  with  the 
development  of  inertial  instruments,  gradiometers,  which  measure 
elements  of  the  matrix.  The  following  discussion  covers  certain 
interesting  mathematical  properties  and  addresses  the  problem  of 
cjqjressing  the  tensor  in  arbitrary  coordinate  frame 

The  representative  gravitation  function,  G^(jr®),  is  only  defined 
for  an  Earth-relative  argument  which  we  indicate  by  the  generic  e-frame 
notation.  Now,  define 


^e(£®)  = d G(r^)  /d  r® 


(A-1) 


This  quantity  will  be  referred  to  as  the  "e-frame  gravity  tensor.  " 

To  explore  some  of  the  mathematical  properties  of  this 
e-frame  gravity  tensor  consider  its  origined.  form.  The  tensor  is  the 
spatial  derivative  of  the  gravitation  function  which  is  in  turn  the  spatial 
derivative,  or  gradient,  of  the  gravitational  potential  function,  V(£®). 
So, 


G®(r«)  = dV(r®)/d  r®  ' (A-2) 

and 
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Te(r®)  = d^V(r®) /(d  r®)^.  (A-3) 

Consider  now  a Cartesian  e-frame  with  orthorgonal  axes  and 


Then, 


(A-4) 


(A-5) 


where  the  subscript  on  V indicates  the  appropriate  partial  derivative. 
We  shall  now  assume  that  the  mass  of  the  air  above  the  Earth's  surface 
can  be  neglected.  In  practice,  gravity  data  is  reduced  taking  this 
assumption  into  account.  Then,  in  the  region  of  interest  V(£®)  is 
harmonic  (Ref  14;  126- 145).  This  property  allows  us  to  state  that  V(£®) 
and  are  continuous.  Then,  the  second  partial  derivatives  com- 

mute, For  example. 


Vjjy  = a^V/axSy  = a^V/dy^x  = 


Since, 


(A-6) 


(A-7) 


From  (A-6)  we  see  that  (A-7)  is  symmetric.  Again,  since  the  potential 
function  is  harmonic,  it  satisfies  LaPlace's  equation: 

= V„  + Vyy  + = 0.  (A-8) 
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Trace 


|re(re)|  = 0. 


(A-9) 


With  symmetry  and  (A-8)  we  conclude  that  of  the  nine  elements  of  the 
gravity  tensor  only  five  are  independent. 

These  mathematical  properties  of  r^£®)  are  tied  to  the 
e-frame.  Since  we  encounter  other  frames  of  reference  in  measure- 
ments and  calculations,  we  need  a means  of  expressing  the  gravity 
tensor  in  these  frames.  This  problem  is  similar  to  transforming  an 
inertia  tensor  from  one  frame  to  another,  and,  not  unsurprisingly,  the 
answer  is  a transformation  of  the  same  form.  Consider  some  arbitrary 
orthogonal  k-frame.  We  can  express 


G^(r®)  = G®(r®)  = cj  G®(C?rk)  (A-10) 


d G^r®) /d  r^  = [d  C®(r®) /d  r®]  [d r®/d r^]  (A-1 1) 


And  we  get. 


r^(r«)  = r^(r')  C J 


(A-12) 


With  these  derivations  and  explanations  the  use  of  the  gravity 
tensor  should  be  clearer. 
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Appendix  B 


Fundamental  INS  Properties 

The  propagation  of  INS  errors  is  modeled  by  first-order 
perturbation  techniques.  References  1 and  2 present  a detailed  develop- 
ment along  this  line.  Some  of  the  more  basic  modes  of  inertial  naviga- 
tion can  be  discerned  from  even  simpler  models.  The  vertical  channel 
of  an  inertial  system  is  unstable,  and  the  horizontal  channels  display  an 
oscillation  known  as  the  Schuler  mode.  These  fundamental  character- 
istics appear,  in  some  form,  in  all  INSs.  The  following  development 
shows  these  underlying  relationships  between  the  position  estimates  and 
the  gravity  model  in  extremely  simple  terms. 

We  have  from  previous  developments: 

+ r (£®)  6r  (9) 

where 

r(r«)  = aC(r®)/ar 

For  convenience,  suppose  at  t " 0 the  Cartesian  z-axis  goes  through  the 
true  position.  The  gravitation  function  can  be  viewed  as  an  expansion 
about  the  centroidal  point  mass  function, 

S£)  = --^  r.  (B-1) 

r^ 
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The  remainder  of  the  gravity  field  will  be  included  in  6£.(£®)  cind  will  be 
ignored  for  this  development  in  order  to  concentrate  on  the  6£  effects. 
We  czm  now  write  much  simplified  error  equations  for  the  vertical  and 
horizontal  channels. 

For  the  vertical  channel,  we  get  from  (9): 


5z  = Sg^  + [3G^/J5K]p6x  + [^Gj^/^y]p  6y  + [^G^/^zjp  6z  (B-2) 

where  point  P for  this  case  will  be  (0, 0,  r^)  for  (x,  y,  z).  In  this  verti- 
cal channel  case,  we  will  assume  8g^,  6x,  and  6y  are  zero.  From 
(B-1)  we  get 

G^  = - p/z^  (B-3) 

So, 

[aG^/az]^^^^  = 2p/r^  (B-4) 

Identifying  p/rQ  as  g,  equation  (B-2)  becomes 


6z  = (Zg/r^)  6z 


(B-5) 


This  homogeneous  error  differential  equation  indicates  solutions  of  the 
form 

6z  = A e"^  + B e"®*^  (B-6) 


where 


a 


and  A and  B are  constants. 


The  positive  exponential  term  demonstrates  the  underlying  instability 
of  the  vertiCcil  channel.  W e should  note  here  that  altimeter  aiding  of 
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the  vertical  loop  is  the  standard  remedy  for  this  INS  problem. 


Now,  turning  our  attention  to  the  horizontal  channel,  we  may 
write  an  equation  in  x analagous  to  (B-2).  This  time  we  will  assume 
6x  is  the  only  active  term  and  get 


6x  = _ q6x 


(B-7) 


The  term  hC^/^x.  can  be  obtained  from  the  following  vector  diagram  if 
we  assume  8x  is  approximately  equal  to  the  chord  in  both  magnitude 
and  direction. 


£,  r^  are  true,  estimated  jsition. 

Cf  6 are  true,  estimated  gravitation. 

0 is  central  angle  in  the  plane  > 
of  estimation  error. 


‘Earth  Center  of  Mass 


By  using  a first  order  approximation  we  get. 


G = G + [ac/^x]p  6x 


(B-8) 


^ * 

From  the  diagram  we  get,  G = G + £ . 


(B-9) 
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So, 


fi.  ” [SG/^x]p6x. 


(B-IO) 


Using  small  angle  approximation  we  get. 


igi  = IGI  *16  6' 


(B-11) 


and 


r 66  = 6 X. 


(B-12) 


Combining  (B-10),  (B-11),  and  (B-12)  yields 


IGI  • 16x1  I , I 

~ dx)pSxj  which  reduces  to 

(1/r)  |G|  = g/r  = |OG/3x)pl  (B-l3) 

Since  ^is  parallel  and  opposite  in  direction  to  the  Sx  line  segment,  we 
conclude  that(3Gx/ax)p  = -g/r.  (B-14) 

Then,  (B-7)  becomes 

8x  = (-g/r)Sx  (B-15) 

From  this  homogeneous  error  differential  equation,  we  expect  solutions 
of  the  form 


X = A sin  w t + B cos  oj  t (B-16) 

8 8 


an  84  min.  period  and  where  A and  B are  constants. 
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-This  mode  was  first  postulated  in  a gyrocompass  design  application 
(Ref  3:15-16),  however  it  appears  in  all  terrestrial  INS  applications 
(Ref  1).  These  characteristic  oscillations  of  the  estimate,  about 
the  true  position,  result  because  horizontal  position  errors  create 
an  opposite  error  in  gravitational  acceleration. 
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Appendix  ^ 


Anomalous  Gravity  Statistical  Methods 

We  know  the  gravity  function  through  a relatively  sparse  sampling 
using  measurements  subject  to  errors.  Our  sketchy  and  imperfect 
knowledge  of  this  field  is  insufficient  for  msiny  analyses  in  Physical 
Geodesy  (e.g.  Ref  65)  and  Inertial  Navigation  (e.g.  Ref  21).  A statisti- 
cal approach  to  these  problems  can  provide  predictions  or  error  assess- 
ments. We  need  a statistical  description  of  the  deterministic  gravity, 
or  anomalous  gravity,  field  to  perform  these  tasks.  The  anomalous 
field  is  completely  specified  by  boundary  conditions  of  either  anomalous 
potential  or  gravity  magnitude  over  a reference  closed  surface.  The 
statistics  are  also  completely  specified  by  the  statistics  of  either  of 
these  quantities  on  the  reference  surface.  We  require  at  least  second- 
order  statistics  for  such  an  analysis.  The  mean  of  either  of  these  basis 
quantities  should  be  zero  since  the  zeroth  term  of  our  model  totally 
accounted  for  Earth's  mass.  The  following  discussion  summarizes 
some  methods  used  to  approximate  or  simulate  the  required  second- 
order  statistic:  the  covariance  function. 

The  basic  "signal"  of  the  anomalous  field  is  represented  by  either 
the  gravity  anomaly  or  anomalous  potential.  The  quantity  required  in 
an  analysis  may  be  a related  term  such  as  deflection  of  the  vertical  or 
geoidal  height  (Ref  66:67).  Figure  C-1  graphically  portrays  these 
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Figure  C-1.  Anomalous  Gravity  Quantities 
(Ref  14:83) 


quantities;  References  14  and  67  give  detailed  definitions  for  these  and 
other  anomalous  gravity  manifestations.  Some  of  the  interrelationships 


Gravity  anomaly  vector: 


(C-1) 


Gravity  anomaly;  Ag  = Bp  ' Vq 


(C-2) 


_Y  is  used  here  for  the  reference  or  normal  field.  This  field  may 
or  may  not  coincide  with  the  model  field 


Geoid  height  or  undulation:  N,  as  shown. 


Anomalous  potential: 

T = YqN 

(C-3) 

Also  relating  Ag  and  T: 

Ag  = - 3T  _ 2T 
3r  r 

(C-4) 

Radial  gravity  disturbance  is  sometimes  used  with  deflection  angles  as 
a complete  definition  for  the  anomaly  vector.  The  radial  disturbance 
is  simply  6g  = — ^ (C-5) 

Through  these  interrelationships,  the  statistics  of  all  these 
anomalous  gravity  manifestations  are  coupled.  It  has  been  shown  (Ref 
66:94-98)  that  the  individual  covariance  functions  for  all  these  quantities 
can  be  derived  from  either  the  anomalous  potential  or  the  gravity  anomaly 
covariance  function.  The  anomalous  potential  covariance  function  is 
mathematically  easier  to  manipulate,  however  the  gravity  anomaly  covar- 
iance function  empirical  approximations  are  more  directly  accessible 
since  gravity  measures  are  the  prime  data.  Once  we  select  the  basis 
quantity,  say  gravity  einomaly,  we  need  to  clearly  define  the  covariance 
relationships  of  these  quantities. 

A covariance  function  must,  in  general,  be  defined  for  a two- 
position-vector  argument.  This  geiierality  should  be  kept  in  mind,  but 
we  do  not  have  data  to  empirically  calculate  such  a function.  We  do  have 
data  on  the  Earth's  surface  and  we  can  estimate  this  function  over  this 
closed  surface.  The  problem  is  simplified  by  reducing  these  Earth 
surface  measurements  to  a surface  which  is  more  mathematically 


convenient:  the  Bjerhammar  sphere  (Ref  68:25-26)  which  approximates 
the  reference  ellipsoid.  This  function  can  be  continued  upward,  if 
necessary,  to  yield  the  function  at  another  level  or  to  show  the  correla- 
tion between  anomalies  at  two  different  levels.  The  upward  continua- 
tion follows  from  the  harmonic  nature  of  the  anomalous  potential.  With 
this  general  background,  we  shall  discuss  more  fully  the  reference 
surface  covariance  function  for  gravity  anomaly  before  introducing  the 
approximation  techniques. 

The  definition  of  gravity  anomaly  covariance  function  is  couched 
in  deterministic  (mean)  rather  than  stochastic  (expectation)  terms. 

Since  the  mean  of  the  gravity  anomaly  over  the  Bjerhammar  sphere  is 
zero  (Ref  14:252),  the  covariance  function  is  the  autocorrelation  func- 
tion (ACF).  By  restricting  to  the  reference  sphere,  we  can  specify  any 
two  points  with  four  arguments.  That  is  C(£j,£2)  can  become 
^2**^2^'  where  X.  is  the  longitude  and  <|>  is  the  latitude.  This  function  can 
be  approximated  for  those  points  with  gravity  measurements,  but  a more 
representative  form  is  needed  to  approximate  the  relationship  between 
unmeasured  points.  Two  forms  of  these  average  or  representative  func- 
tions appear  in  the  literature  most  often. 

The  first  form,  C(i)^,«),  is  defined  in  terms  of  a shift  (arc) 
distance  s = R4<  and  an  arc  azimuth  argument  o . Here,  ^ represents 
the  central  angle  between  the  two  points  on  the  sphere.  Figure  C-2 
portrays  this  relationship.  The  defining  equation  is 

C(4'» «)  = [Agp  Agpi] 
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(C-6) 


where  is  the  mean  over  the  set  K which  is  defined  by 
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WiT  I vtn  I 
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(1)  P and  P'  are  on  the  sphere; 


K = 


(P,P‘) 


(2)  Central  angle 
and 


is  ijj ; 


(3)  Arc  PP  has  azimuth  a . 


C(4<, a)  is  then  the  average,  or  mean,  pair  product  over  all  point  pairs 
on  the  sphere  which  subtend  a central  angle  of  ij*  and  which  have  a join- 
ing arc  of  azimuth  a.  This  average  covariance  function  is  consistent 
with  the  assumption  of  homogeneity  (Refs  66:85  and  69:1). 

The  other  common  form  of  the  covariance  function  is  generated 
in  a manner  that  is  consistent  with  symmetry  in  the  argument  4^.  This 
radial  (in  the  horizontal  distance  sense)  symmetry  is  a characteristic 
of  an  isotropic  field  over  the  sphere  (Refs  66:85  and  69:1).  The  function 
C(i|J)  is  easily  defined  in  terms  of  0(4*,“  ) by 


C(4^)=  M [C(4>,«)]. 

ttc(0,2ff) 


(C-7) 


The  assumption  of  an  isotropic  field  may  not  be  valid  for  a particular 
section  (see  Ref  49:37-39  and  Ref  70:26  for  example  real  data),  but  this 
assumption  is  attractive  for  its  mathematical  simplifications  (Ref  71:39). 

Equation  (C-7)  can  be  approximated  with  available  gravity  measure- 
ment data  which  results  in  a tabular  form  of  C(ip),  Examples  of  such 
tabulations  abovmd--each  generated  under  different  rules  or  data  base 
(e.  g.s  Ref  14:254  and  Ref  69:80-83).  Such  a table  is  our  only  true 
measure  of  the  covariance  function,  however  this  form  is  cumbersome 
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for  mathematical  manipulations.  A closed  form  mathematical  approxi- 
mation is  useful  in  many  applications,  and  many  such  approximations 
have  been  suggested  (Ref  72:10-13  and  Ref  73:10-14).  The  Hirvonen 
covariance  function  model  (Ref  14:255)  is  one  example 

C 

where  R is  spherical  radius  and  and  d are  parameters.  In  the  cited 
example,  a value  of  337  mgal  was  provided  for  Cq,  and  d was  speci- 
fied as  40  km.  These  parameters  were  selected  for  a specific  geo- 
graphic area  (southern  Ohio). 

The  selection  of  such  model  parameters  is  critical  for  the  uses 
of  these  closed  form  approximation.  These  expressions  are  primarily 
used  for  statistical  weightings  in  the  region  of  small  shift  distance 
(within  the  half -variance  range,  see  Figure  C-3),  so  the  behavior  of 
these  closed  form  approximations  in  this  region  receives  special 

Co 


Figure  C-3.  Covariance  Shaping  Factors 
(Ref  72:22) 
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attention.  The  three  shaping  factors  which  are  used  to  judge  the  good- 
ness of  fit  of  these  approximations  are  (1)  the  variance,  C^,  (2)  the 
radius  of  curvature,  p,  at  zero  shift,  and  (3)  a representative  distance, 
d,  quantity  defined  by  the  shift  distance  to  the  half-variance  level. 

Figure  C-3  portrays  these  parameters;  Reference  72  has  a thorough 
discussion  of  these  issues  in  Section  3. 

Another  modeling  approach  is  to  develop  a closed  form  approxi- 
mation for  the  spherical  harmonic  coefficients  of  the  covariance  func- 
tion spherical  harmonic  decomposition.  These  coefficients  are  referred 
to  as  degree  variances  since  they  represent  the  variance  of  the  particu- 
lar wavelength  components.  This  technique  is  developed  in  Reference 
69.  This  xnodeling  technique  is  sometimes  combined  with  the  previous 
technique:  the  closed  form  model  is  used  for  the  near  zero  shift  region 
and  the  degree  variance  model  for  the  remainder  of  required  covariance 
evaluations . 

Our  discussion,  thus  far,  has  been  confined  to  the  reference 
surface  form  of  the  covariance  function.  The  more  general  form  which 
allows  points  P and  P'  to  be  at  different  levv'jls  needed  for  some 
analyses.  An  example  of  this  need  to  span  altitude  differences  is  in  the 
least-squares  combination  of  satellite  data  with  surface  data  (Ref  74). 

If  we  construct  the  covariance  function  with  arguments  and  we 
can  theoretically  continue  the  function  upward  in  one  argument,  say  ££» 
while  fixing  the  other,  to  the  reference  surface  (e.  g.  Ref  68:32). 
Heller  and  Thomas  (Ref  50:Appendix  B)  have  recently  demonstrated 
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some  concise  analytical  techniques  to  haindle  such  upward  continuations 
in  a more  general,  spherical  rather  than  planar,  manner. 

These,  then,  are  some  of  the  statistical  modeling  approaches  in 
Physical  Geodesy.  Inertial  Navigation  analyses  are  interesting  in  com- 
parison. The  Physical  Geodesy  task  is  based  on  sample  measurements 
of  anomalous  gravity  and  a statistical  model  used  in  conjunction  to 
adjust  data,  predict  gravity,  or  estimate  the  associated  errors.  The 
inertial  analyses  attempt  to  characterize  INS  errors  based  on  a statisti- 
cal model  of  gravity  modeling  errors  and  a dynamic  error  propagation 
model.  The  tasks  while  quite  different  share  the  covariance  statistics. 

A covariance  function  or  autocorrelation  function  (ACF)  for 
gravity  anomaly  io  useful  in  these  INS  statistical  error  analyses.  An 
INS  is  expected  to  operate  globally — or  at  least  anywhere  within  some 
large  region.  We  wish  to  characterize  INS  performance  throughout  this 
region.  If  the  cuiomalous  gravity  field  were  known  exactly  throughout 
this  area,  we  would  need  data  from  the  ensemble  of  possible  trajectories 
to  calculate  average  INS  error  statistics.  Fortunately,  given  an  ACF,  • 
we  can  estimate  these  INS  performance  indices  by  simulating  the 
anomalous  field  with  stochastic  processes.  These  processes  can  be 
cast  in  a form  compatible  with  the  linear  time -invariant  INS  error 
propagation  models  (Refs  1 and  2).  The  subject  ACF  is  the  key  to  this 
modeling  process. 

The  deterministic  anomalous  gravity  field  for  some  area  is 
characterized  by  second  order  statistics--means  and  covariances --of 
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the  gravity  anomaly.  The  INS  error  performance  is,  then,  simulated 
by  the  performance  of  an  INS  driven  by  the  outputs  of  a shaping  filter 
driven  in  turn  by  white  Gaussian  noises  (Ref  75:217-218).  The  ensemble 
of  trajectories  is  simulated  by  the  ensemble  of  noise  functions  in  that 
sample  space.  The  shaping  filter  design  and  the  strength  of  the  associ- 
ated noises  are  selected  such  that  the  filter  output  covariance  approxi- 
mates well  the  empirically  determined  (estimated)  anomaly  ACF.  The 
mission  velocity,  as  we  have  previously  seen,  translates  these  spatial 
ACFs  to  temporal  ACFs. 

This  concept  of  approximation  has  several  years  of  development 
history.  The  conflicting  requirements  of  mathematical  tractability  and 
physical  reasonableness  have  channeled  this  effort.  We  shall  identify 
some  major  areas  of  interest  in  this  evolutionary  process  before  pro- 
ceeding with  a chronological  account  of  the  milestones  in  the  develop- 
ment of  this  modeling  concept. 

Most  gravity  data  is  reduced  to  either  the  reference  ellipsoid  or 
the  geoid,  so  we  can  only  readily  calculate  an  ACF  for  trajectories 
bound  to  these  surfaces.  These  data  are  in  the  form  of  gravity  anomaly 
or  geoid  vindulation  (anomalous  potential  is  linearly  related  to  geoidal 
undulation  by  Brun's  formula',.  The  practical  problem  is  to  statistically 
produce  the  INS  error  inputs  at  mission  altitude  and  velocity  based  on 
these  available  statistics.  The  physical  interrelationships  between 
these  simulated  inputs --vertical  deflections,  gravity  disturbance,  and 
geoid  undulation- -must  be  consistently  reflected  in  the  joint  statistics 


of  these  inputs.  A concise  method  £or  upward  continuing  these  data  to 
mission  altitude  is  desired,  and  as  mentioned,  the  model  must  blend 
with  INS  models  for  maximum  usefulness. 

The  fundamental  modeling  choices  are  (1)  which  anomalous  grav- 
ity variable  to  use  as  a basis  for  deriving  other  ACFs  and  (2)  which  var- 
iables to  use  for  INS  error  inputs.  Early  studies  were  quasi-static 
parametric  analyses  of  the  horizontal  channels  of  a local-level  INS 
mechanization.  The  vertical  deflection  angles  were  natural  choices 
for  this  task.  Later  studies  generalized  the  deflection  angles  to  “along- 
track"  and  "cross -track"  coordinates  to  fully  exploit  tlie  isotrooic 
assumption.  Barometric  altimeter  aiding  of  the  vertical  channel  intro- 
duces an  additional  disturbance  error  driving  source  due  to  geoidal 
undulation.  Since  all  these  sources  can  be  related  through  surface 
integrals  to  either  anomaly  or  undulation,  it  is  physically  unreasonable 
to  assume  them  all  to  be  statistically  uncori elated.  Answers  to  these 
modeling  questions  have  evolved  over  the  years.  A chronological 
account  of  this  development  is  a suitable  context  for  further  discussion. 

Levine  and  Gelb  (Ref  21)  developed  models  to  approximate  ACFs 
of  the  prime  and  meridional  vertical  deflections.  The  physical 
processes  were  simulated  as  independent  outputs  of  first-order  spatial 
stochastic  differ.ential  equations  of  the  form 

I W+q^(x)  (16) 


120 


where  x is  shift  distance,  is  correlation  distance,  and  is  a white 

2 2 

Gaussian  noise  process  with  power  spectral  density  of  2cr|/d.  is 
then  the  process  output,  hence  vertical  deflection,  variance.  Tne 
parameters  <r^  and  d^  (Ref  76:86-88)  can  be  varied  to  produce  the  best 
overall  fit  between  the  empirically  derived  ACF  and  the  exponential  ACF 
resulting  from  (16): 


(14) 


This  simple  model  was  criticized  on  two  counts.  First,  the  ACF  is 
non-negative  and  empirically  derived  ACFs  have  evinced  negative 
regions.  We  shall  see  later  that  a positive-only  covariance  function  is 
physically  unreasonable;  in  fact  it  must  be  zero  mean  over  the  sphere. 
Also,  the  ACF  is  expected  to  be  rounded  (zero  slope)  at  zero  shift.  The 
surface  integral  relationships  of  anomaly  to  Earth  mass  distribution 
implies  at  least  first-derivative  continuity  at  zero  shift.  This  notion  is 
also  observed  in  empirically  generated  anomaly  ACFs,  but  is  lacking  in 
the  first-order,  exponential  model. 

These  physical  reasonableness  objections  can  be  overcome  at  the 
cost  of  additional  model  mathematical  complexity.  Shaw  et  al  (Ref  77) 
led  the  way  by  forming  more  complete  models.  They  developed  the 
relationships  between  deflection  ACFs  which  we  want  and  anomaly  ACF 
which  we  have- -at  least  empirically.  Since  veitical  deflections  can  be 
calculated  from  anomaly  data  by  the  Vening  Meinesz  formula  (Ref  14: 
114),  this  relationship  provides  the  link  between  anomaly  ACF  and 
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deflection  ACFs  and  deflection  cross -correlation  function  (CCF).  Flat 
Earth  approximations  are  used  along  with  the  assumptions  that  the 
gravity  anomaly  is  homogeneous  and  isotropic  on  the  geoid.  The  anom- 
aly ACF  can  then  be  expressed  in  terms  of  a horizontal  radial  shift 
distance.  The  Bessel  function  which  follows  quite  naturally  is  not  com- 
patible with  our  linear  filter  goal.  The  exponential  ACF  of  Reference 
21  (Levine  and  Gelb)  is  further  developed  in  Reference  77  by  the  intro- 
duction of  deflection  CCF.  An  improvement  was  identified  in  the  CCF. 

It  was  shown  chat  a change  of  coordinates  from  the  traditional  east-west 
and  north-south  for  deflection  definition  could  result  in  a zero  CCF.  If 
generic  "along-track"  and  "cross-track"  coordinates  are  used  for 
expressing  deflection  components,  the  isotropic  assumption  then  implies 
zero  deflection  CCF.  Since  the  partial  derivatives  used  in  defining  the 
deflection  angles  are  not  tied  to  a particular  azimuth,  the  statistical 
representation  is  an  average  over  all  azimuths  which  yields  this  desir- 
able modeling  result. 

Kasper  (Ref  78)  added  another  3tep  with  second-order  Gauss- 
Markov  process  models.  He  achieved  the  zero-slope-at-zero-shift 
goal,  but  continued  v/ith  a positive-only  ACF.  Simulation  results  com- 
paring this  second-order  model  to  the  previous  first-order  one  show 
almost  identical  rms  position  error  per  unit  of  input  standard  deviation. 
The  rms  velocity  errors  were  apart  by  a factor  of  two,  however.  The 
lower  errors  associated  with  the  second-order  Markov  covariance 
model  compare  favorably  with  results  obtained  with  measured  gravity 


122 


disturbances  as  INS  input.  This  result  demonstrates  the  importance  of 
the  rounded  correlation  feature.  This  should  come  as  no  surprise  since 
the  importamce  of  the  near-zero-shift  region  received  so  much  attention 
in  the  discussion  of  Physical  Geodesy  applications. 

Bellaire  (Ref  79)  developed  methods  for  continuing  these  geoid  or 
ellipsoid  based  ACFs  up  to  mission  altitude.  By  employing  harmonic 
function  analysis,  the  upward  continuation  is  cast  in  the  form  of  surface 
integrals.  The  anomaly  ACF  is  shown  to  be  non-stationary  in  altitude 
due  to  attenuation.  The  anomaly  field  is  stationary  in  the  horizontal 
arguments  at  all  levels  if  it  is  so  on  the  geoid.  The  ACF  between 
anomalies  at  two  different  altitudes  is  shown  to  be  equal  to  the  ACF  at 
the  average  altitude.  This  new  analytical  tool  allows  a much  more 
general  trajectory  study  than  the  constant-altitude,  constant-velocity 
previous  approach. 

The  models  and  methods  of  Kasper  (Ref  78)  and  Bellaire  (Ref  79) 
support  a statistically  coordinated  simulation  of  anomalous  gravity 
driving  all  three  INS  axes  over  a range  of  variable  altitude  and  velocity 
trajectories.  The  only  important  unmodeled  effect  at  this  point  is  the 
geoidal  undulation  disturbance  of  the  barometric  altimeter  used  to  aid 
the  vertical  channel.  Jordan  (Ref  80)  provided  the  theoretical  develop- 
ment to  include  this  term,  ivlor-'  importantly,  however,  Rafeience  80 
introduces  a set  of  mathematical  consistency  "necessary  conditions" 
with  which  we  can  judge  these  proposed  models. 
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These  consistency  criteria  are  based  on  our  understanding  of 


anomalous  potential  and  its  relationships  to  gravity  anomaly  and  to 
vertical  deflections.  The  method  consists  of  starting  with  the  ACF 
which  our  statistical  simulation  produces,  say  the  anomaly  ACF.  The 
undulation,  or  anomalous  potential,  ACF  is  derived  from  this  anomaly 
ACF  using  the  relationship 

(iT2 

where  is  undulation  ACF,  <j>  is  anomaly  ACF,  x and  y are  hori- 

oo 

zontal  plane  shift  distances,  and  Vq  is  reference  gravity.  We  can  relate 
the  second  partials  above  to  the  negatives  of  the  vertical  deflection  ACFs. 
Since  these  deflection  ACFs  at  zero  shift  represent  variances  they  must 
be  positive  there.  Hence,  a derived  undulation  ACF  which  is  not  con- 
cave at  zero  shift  indicates  the  original  anomaly  ACF  is  physically 
unreasonable.  Since  the  kernel  associated  with  (C-9)  is  singular, 
another  specific  condition  must  be  met  by  the  anomaly  ACF  for  the 
associated  undulation  ACF  to  be  bounded  everywhere  on  the  plane.  The 
constraint  is  that  the  anomaly  ACF  must  be  zero  mean  over  the  plane. 

In  summary,  if  our  statistical  model  generates  an  anomaly  ACF  which 
yields  an  associated  undulation  ACF  which  is  either  unbounded  or  is  not 
concave  down  at  the  origin,  that  statistical  model  is  physically  incon- 
sistent. 

The  previous  models  were  each  found  lacking  in  some  way.  New 


second-  and  third-order  Gauss-Markov  process  models  were  proposed. 
Even  this  second-order  model  fails  to  meet  all  criteria.  Both  new 
models  accommodate  a negative  ACF  at  some  shift  distance,  but  only 
the  third-order  model  has  the  rounded  zero-shift  profile.  This  third- 
order  model  represents  the  apex  of  this  evolution. 

What  more  could  we  ask  for?  The  models  discussed  above  all 
assume  a flat  Earth  using  the  justification  that  correlation  distance  is 
much  smaller  than  the  Earth's  radius.  Statistical  studies  of  long-range 
missions,  such  as  ICBM  trajectories,  require  a more  global  setting. 
Heller  and  Thomas  (Ref  50)  have  begun  v/ork  on  such  modeling  tech- 
niques. They  propose  modeling  the  globai  free-air  anomaly  field  as 
the  result  of  buried  spherical  white  noise  processes  which  are  continued 
up  to  and  summed  on  the  reference  surface.  This  technique  has  success- 
fully modeled  some  empirical  anomaly  ACFs.  The  present  approach, 
however,  yields  an  all-positive  ACF.  So,  research  on  the  global  model- 
ing problem  is  likely  to  continue. 

In  summary,  we  have  seen  the  development  and  use  of  the  gravity 
anomaly  ACF  or  covariance  function  from  two  different  points  of  view. 
The  Physical  Geodesy  problems  require  the  covariance  function  as  an 
interpolation  or  prediction  tool.  The  Inertial  Navigation  studies  require 

The  term  "Third-order  Markov"  means  that  either  the  anomaly 
or  undulation  can  be  modeled  as  a Third  order  linear  system  driven  by 
White  Gaussian  noise.  To  form  a complete  system  to  statistically  pro- 
duce T),  N,  andAg  requires  eight  integration  levels  (Ref  81:3-29  to 
3-31). 
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the  anomaly  ACF  to  supply  the  input  parameters  to  stochastic  models 
which  allow  simulation  of  INS  errors  we  expect  the  anomalous  field  to 
produce. 
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Appendix  D 


INS  Gravitational  Model  Performance 
Cost  Function 


Gravitational  modeling  errors  should  be  studied  in  terms  of 
their  effect  on  mission  performance.  For  example,  we  might  require 
the  overall  system  to  meet  some  circular-error-probable  (CEP)  cri- 
teria. The  performance  cost  of  each  system  component  would  then  be 
stated  in  terms  of  the  effect  on  CEP.  In  such  statistics,  it  is  custom- 
ary--and  justifiable  in  most  cases--to  consider  the  overall  system 
performance  estimate  to  be  the  root- sum-square  (rss)  of  the  individual 
(supposedly  independent)  contributions.  The  effects  of  gravitational 
modeling  errors  can  be  statistically  characterized  in  terms  of  the  INS 
error  state  covariance  matrix,  P„„(t).  This  error  covariance  matrix 
can  be  processed  to  yield  the  CEP  cost  due  to  the  gravitational  model- 
ing errors  alone.  With  this  example  as  justification,  the  following 
presentation  describes  how  can  be  processed  to  form  such  a 

system  cost. 

We  begin  by  discussing  a terminal  cost  strategem.  In  many 
instances  of  payload  or  weapon  release,  the  t ;tal  delivery  accuracy 
can  be  related  to  the  terminal  state  estimation  accuracy  by  a linear 
propagation  of  errors  using  linear  perturbation  techniques.  Let  ^ be 
the  system  "miss"  coordinates;  let  X(t)  be  the  INS  error  state  vector 
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{e.  g.  downrange  and  crossrange);  and  assume  terminal  error  is  given 


by 

l = IK(tf)  (D-1) 

where  tj  is  terminal  time. 

We  can  form  a terminal  cost  by 

J = y'^y  = Trace  I yy*^!.  (D-2) 

To  be  more  general  we  can  weigh  the  "miss"  vector  with  a matrix  W 
which  can  be  considered  symmetric  without  loss  of  generality.  Such 
a W matrix  should  be  either  positive  semi-definite  or  positive  definite. 
For  terminal  cost  we  get 

J = Z^(t^)Wz(tf)  = xT(tf)H%HX(t^) 

= Trace  jx(t^)x’^(t^)H'^WHl  (D-3) 

Since  our  error  state  X(*)  is  being  treated  statistically,  we  are  inter- 
ested in  the  expectation.  S',  of  J over  the  ensembla  in  our  sample 
space: 

S[J]  =(?)  Trace  [X(tf)X’^(tj)H'^WH]  | 

= Trace  j^[X(tf)x'^(tf)]H’^WH} 

= Trace  (D-4) 
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assuming  H and  W are  invariant  over  our  sample  space.  An  analagous 
running  cost  function  can  be  defined  for  those  missions  which  require 
INS  accuracy  over  all  mission  time.  Again,  we  can  conceive  of  a 
transformation  to  those  coordinates  on  which  system  performance  is 
judged. 


i:(t)  = H(t)X(t) 


(D-5) 


Then 


j(t)  = l'^(t)W(t)x{t) 

= Trace  j X{t)x'^(t)  H'^(t)W(t)H(t)(  (D-6) 


Then, 

J(t)  = y'Trace|X(p)x'^(p)HT(p)W(p)H(p)|dp  (D-7) 

o 

And 

^[J(t)]  = f Trace|p^^(p)H'’^(p)W(p)H(p)|dp  (D-8) 

0 


assuming  H(t)  and  'W(t)  are  invariant  over  our  sample  space.  The  cost 
functions  given  by  (D-4)  and  (D-8)  are  reminiscent  of  stochastic  termi- 
nal controller  and  stochastic  regulator  cost  functions  (Ref  82:415). 
Indeed  the  model  parameters  involved  in  the  optimization  process  may 
be  viewed  as  control  variables  we  use  to  minimize  these  costs.  As 
with  stochastic  control  problems,  there  may  be  missions  where  a 
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combination  of  terminal  and  running  costs  are  appropriate.  The 
fundamental  quantity  needed  for  these  cost  evaluations  is  the  INS 
error  state  covariance  matrix,  ^xx^  t)  • 
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Appendix  E 


The  Q- Matrix 


The  Q-matrix  statistically  summarizes  the  anomalous  gravita- 
tional information.  As  mentioned  in  Appendix  C,  the  covariance  func- 
tions of  various  anomalous  gravitation  quantities  can  be  derived  from 
either  the  anomaly  covariance,  C,  or  the  anomalous  potential  covari- 
ance, K.  This  result  holds,  as  well,  for  the  residual  field  after 
inclusion  of  the  disturbance  model  Sgm  since  this  model  is  based  on  a 
harmonic  potential  field.  So,  while  the  following  discussion  uses  the 
gravitation  disturbance,  Sg,  and  the  anomalous  potential,  T,  the  same 
relationships  hold  for  residual  disturbance,  and  residual  potential, 
8T.  By  our  previous  definition  we  have 


Q[r.,r2;K(R,R«;0)]  = S 


8£(li)S£.'^(l2) 

T{ri)S£T(^) 


S£.(£l)T(r2)' 
T(ri)T(r2)  _ 


(E-1) 


where  and  £2  ^ design  mission  trajectory  Let  x, 

y and  z be  our  generic  Cartesian  coordinates;  then. 


Sgx(£l) 

8gy(r  1 ) 

8gz(£i) 


(E-2) 


and  Sgjj(rj)  = 


J.1. 

3xi 


(E-3) 
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With  the  understanding  from  (E-3)  notation  (E-1)  can  be  rewritten  as 


Q[r,,r2;K(R.R';e)]  = (^ 


* T(rj)T(r2) 


’ 82 
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(E-4) 


For  practical  reasons,  we  expect  the  partial  deriv^^tives  to  be  uni- 
formly bounded.  Hence,  these  operations  are  bounded  and  linear. 

So,  we  may  interchange  the  partial  derivatives  with  the  linear  expecta- 
tion operator.  Then,  (E-4)  becomes 


Qlr^.r^;  K(R,R';0)]  = 


• a2 

a2 

a2 

a 

axidx2 

ax^ay2 

axj^^zj 
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(E-5) 


Since  K(£j^,r^2l0)  can  be  derived  from  K(^, ^';0),  this  equation  demon- 
strates that  all  necessary  statistics  are  embodied  in  the  single  reference 
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surface  covariance  function.  With  this  point  made,  we  shall  not 
implement  (E-5)  directly.  Rather,  the  elements  of  the  Q-matrix  can 
be  derived  (Ref  50)  from  the  (E-5)  relationship  with  proper  coordi- 
nate choice. 

The  usual  form  of  the  covariance  functions  (Appendix  C)  is  in 
terms  of  a central  angle,  argument  That  is, 

K(R,R>)  = K(i}f)  (E-6) 

where 

^ = cos"^  [R-2  R-R']  . (E-7) 

We  have  dropped  the  sample  space  argument,  here,  since  the  develop- 
ment is  generally  applicable.  It  has  been  shown  (Ref  50)  that  potential 
covariances  of  the  form  (E-6)  can  be  continued  upward  with  arguments 
rj,  r2  and  ip.  That  is, 

^(£^,£2)  = K(rj,r2,(i»)  (E-8) 

A particularly  useful  coordinate  choice  in  this  case  is  Inplane- 
Transverse-Radial  (ITR)  coordinates  shown  in  Figure  E-1. 

If  we  use  the  autocorrelation  notation  of  Appendix  C we  can 
summarize  (E-5)  as 
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Q[tj,,  r2^,i}r,  K(V;d)  \ - 


) <Pxz^  ^ ^ 

^zy^  ^ ^zz^  ^ ^zT^  ^ 

^Ty(  ^ ^Tz^  ^ ^TT^  ^ 


(E-9) 


The  tasks  remaining  are  to  provide  the  various  auto-and  cross -corre 
lation  functions  based  on  the  reference  surface  version  of  the  lower 


right-hand  element  of  (E-9):  R» 


Note:  View  is  in  the 


£2  plane 
z-axis  radial 


x-axis  inplane  (downrange) 


y-axis  transverse 
(crossrange) 


Figure  E-1.  Inplane-Transverse-Radial  Coordinates 
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